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CHAPTEB  i 

INTHO.0VJ'.;  ' ..".'..  . 

During  the  past  decade  and  &   naif,  process  dynamics  and 
control,  specifically  optimal  Control,  has  been  a  subject  of 
considerable  interest  to  engineers  in  all  fields.   Basically  an 
optimal  control  problem  can  be  stated  as  follows i  given 
equations  which  describe  a  dynamic  system,  determine  how  the 
Systems  controls  must  be  manipulated  as  a  function  of  time  in 
order  to  maximize  or  minimize  a  performance,  criterion  subject 
to  constraints  on  the  controls  and/or  on  the  elements  of  the 
system  itself. 

The  problem  of  optimal  control  of  distillation  columns  has 
been  of  continued  interest  to  chemical  engineers  in  recent  years 
and  all  sorts  of  sophistication  from  the  most  recent  developments 
in  automatic  control  theory  has  been  employed  in  solving  these 
problems.  However,  nearly  all  of  these  investigations  have  been 
based  on  ideal  models  with  the  assutnptiorj  of  complete  mixing  in 
the  liquid  phase  on  each  tray.   And  this  In  spite  of  the 
development  of  several  models  for  the  representation  of  liquid 
phase  mixing  on  distillation  tray.?  (see  the  book  by  Holland  [l]). 
The  present  -Study  was  motivated  by  the  belief  "tHat  the  effect 
of  liquid  phase  mixing  on  the  optimal  control  of  distillation 
columns  has  not  been  sufficiently  highlighted  in  the  past.   This 
thesis  is  intended  to  be  a  contribution  toward  obtaining  a  better 
understanding  of. the  problems  involved  in  obtaining  an  optimal 
control  policy  for  a  distillation  column  wherein  the  liquid  phase 


mixing  on  each  tray  Is  described  t:j  the  so-called  mixing  pool 
model  rather  than  the  conventional  completely  mixed  tray. 

The  system  studied  in  this  work  is  a  two  tray  column 
separating  a  binary  mixture  wherein  the  stlllpot  acts  ss  the 
bottom  tray  while  the  top  tray  serves  as  a  rectifier.  The 
special  feature  of  this  column  is  that  the  top  tray  is  not 
considered  as  an  ideal  tray  with  complete  mixing,  instead  it  is 
considered  as  being  divided  into  stages  or  pools,  each  pool 
consisting  of  completely  mixed  liquid.   This  is  Kirschbaura's  [2] 
assumption  that  there  are  pools  of  liquid  on  the  tray  which  are 
completely  mixed.  The  so-called  mixing  pool  model  is  described 
as  being  divided  into  a  number  of  completely  mixed  stages  or 
pools  [3,  *0«   The  concentration  gradient  across  the  tray  is 
assumed  to  be  made  up  of  a  series  of  pools  with  each  mixing  pool 
having  a  uniform  composition.  For  a  very  large  number  of  very 
small  pools,  the  concentration  gradient  becomes  continuous  and 
plug  flow  exists.   Conversely  for  a  tray  that  Is  completely  mixed, 
there  is  no  concentration  gradient  and  the  entire  tray  has  a 
uniform  composition.  These  two  limiting  cases,  however,  may  not 
be  realized  in  practice.   The  effect  of  the  degree  of  mixing 
indicated  by  the  number  of  mixing  pools,  on  the  control  polipy  of 
the  column  -le->the  chief  purpose  of  this  investigation. 

For  the  sake  of  comparison,  a  second  system  which  Is 
referred  to  as  the  "reference  system"  Is  also  considered  in  this 
work.   It  is   identical  to  the  system  under  consideration  except 
for  the  top  tray  which  is  a  conventional  ideal  tray  with  complete 
mixing,  i.e.  with  r.  single  pool.   This  system  serves  as  a  basis 


for  comparison  In  this  Investigation. 

The  procedure  used  here  Is  to  consider  a  single  specific 

I. 
example  wherein  the  system  suffers  a  disturbance  through  the 

feed  composition  which  in  turn  displaces  the  overhead  distillate 

composition  from  its  steady  state,  and  then  to  apply  Pontryagtn's 

Maximum  Principle  in  order  to  obtain  the  control  policy  which 

villi  return  the  overhead  composition  to  its  steady  state  value 

in  the  shortest  possible  time  (time  optimal  problem).  An 

extension  is  also  made  to  the  case  where  the  deviation  from  the 

steady  state  composition  is  kept  at  a  minimum.  We  assume  that 

control  is  effected  by  manipulation  of  the  overhead  reflux  flow 

rate. 

The  basic  theoretical  reference  on  the  Maximum  Principle  is 
thl1  book  by  Pontryagin  and.  his  co-workers  [5].  A  good 
elementary  account  of  the  Maximum  Principle  in  both  continuous 
and  discrete  form,  along  with  numerous  engineering  applications 
can  be  found  in  the  two  works  by  Fan  et  al  [6,  7  J,   While  many 
significant  applications  of  the  Maximum  Principle  have  been 
demonstrated  in  other  engineering  fields,  it  is  only  during  the 
past  few  years  that  chemical  engineers  have  begun  to  apply  this 
mathematical  theory  to  problems  of  chemical  engineering 
importance.  * 

One  of  the  first  and  most  extensive  of  all  studies  In  this 
area  was  accomplished  by  Ciebenthel  and  Aris  [8],  They  treated 
the  optimal  regulation  of  a  continuous  stirred  tank  reactor. 
Coward  [9]  In  his  work  on  the  time  optimal  problem,  is  among 
the  early  investigators  of  the  Maximum  Principle  as  applied  to 


I 

■ 

distillation. 

In  this  study  the  Maximum  Principle  coupled  with  certain 
numerical  techniques  is  employed  for  the  purpose  of  determining 
the  optimal  control  policy.  The  mathematical  models  used  are 
adequate  to  represent  the  phenomenon  being  investigated  without- 
having  to  go  into  the  complexities  of  either  the  hydrodynamics 
or  the  energy  balance  of  the  system. 

Chapter  2  deals  with  the  analysis  of  the  reference  system 
and  here  a  control  policy  is  arrived  at  for  the  regulation  of 
the  overhead  reflux  rate  for  the  time  optimal  problem.   In 
Chapter  3  a  similar  analysis  is  carried  out  with  the  mixing  pool 
model  with  two  pools  or  tanks  in  series,  wherein  the  parameters 
retain  the  same  values  as  those  in  Chapter  2  except  for  the 
introduction  of  the  mixing  pool  model  in  place  of  the  conventional 
conrpletely  mixed  tray.   Once  more  a  control  policy  is  derived  for 
the  regulation  of  the  overhead  reflux  rate  for  the  time  optimal 
problem.  In  Chapter  4  an  extension  is  made  to  the  case  where 
the  control  is  optimum  In  the  sense  of  minimum  total  deviation. 
First  the  control  policy  is  obtained  for  the  reference  system 
and  then  for  the  nixing  pool  model  system.  Chapter  5  summarizes 
the  various  results  from  Chapters  2,  3  and  k   and  these  results 
are  then  analyzed  and  conclusions  drawn.   Also  Several 
supporting  appendices  deal  with  the  analog  computer  flow  charts 
and  digital  computer  programs,  derivation  of  the  Maximum 
Principle  algorithm  and  steady  state  analyses  of  the  systems 
considered  In  Chapters  2,  3  and  t-. 
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CHAPTE8  2 

THE  HEFEHENCE  SYSTEM 

In  this  chapter  the  reference  system  is  considered. 
Initially  the  performance  equations  of  the  system  are  derived. 
The  exact  nature  of  the  problem  stated  and  finally  the 
application  of  the  Maximum  Principle  and  the  derivation  of  the 
optimal  control  policy  terminates  the  chapter. 

1.   DEVELOPMENT  OF  PERFORMANCE  EQUATIONS. 

The  assumptions  and  column  characteristics  involved  in  the 
derivation  of  performance  equations  are  listed  below. 

(a)  A  column  with  two  theoretical  trays  including  the 
Stlllpot  and  using  a  total  condenser,  is  to  be  considered. 

(d)  The  column  is  separating  a  biliary  mixture  assuming  a 
linear  vapor-liquid  equilibrium  relationship,  i.e.,  y  ^  mx,  4   c. 

(c)  The  top  tray  and  the  stlllpot  show  constant  molal 
holdup  which  is  invariant  with  time. 

(d)  There  is  constant  molal  overflow,  i.e.,  liquid  and 
vapor  flows  from  tray  to  tray  are  constant. 

(e)  Vapor  holdup  Is  assumed  negligible. 

(f)  Adequate  cooling  water  and  steam  are  available. 

(g)  To  pel-silt  ample  boiling  surface,  the  stlllpot  holdup 
is  assumed  to  be  2.5  times  the  holdup  on  the  top  tray. 

(h)  The  feed  is  a  completely  condensed  liquid  at  saturation. 

The  essential  features  of  the  system  along  with  the  various 
flow  streams  are  shown  in  Figure  1.  The  feed  is  introduced 


top  troy 


HT  ,  X| 


F,  X, 


condenser     y 


v>y, 


/) 


L,  XD 


L,X, 


H  B  .  X 


B,  A  2 


v,y2 


bottom  troy 


B,X2 


Fig. I.    Distillation    column    with   top  tray   described 

by   a    single-;  completely   mixed    tank  (Reference- 
System)  . 
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directly  Into  the  stillpot  as  a  condensed  saturated  liquid.  In 
the  stillpct  heat  is  transferred  to  the  liquid  and  the  vapor 
leaving  the  stillpot  passes  upward  through  the  rectifying  section 
which  consists  of  a  single  theoretical  tray  (a  single  completely 
mixed  tray).  The  overhead  vapor  is  completely  condensed  and 
part  of  the  resulting  liquid  is  drawn  off  as  overhead  product 
while  the  rest  of  the  liquid  is  returned  to  the  top  tray  as 
reflux.  Liquid  also  leaves  the  reboiler  and  is  drawn  off  as 
bottoms  product. 

On  the  basis  of  the  above  assumptions  and  description  the 
mathematical  relationships  between  the  various  variables  can 
now  be  derived.  The  dynamic  behavior  of  the  column  in  the 
transient  state  may  be  written  in  the  differential  form  by 
performing  a  light  component  balance  on  each  tray  along  with 
overall  material  and  light  component  balances  over  the  entire 
system  as  follows}  (see  Figures  2  and  3). 

For  the  top  tray,  the  light  component  balance  gives 
Input  -  Output  n  Accumulation 

Lxfi  +  Vy2  -  LXl  -  Drx.  HT  -2 


Since 


>■;': 


XD  =  yl  =  mXl  +  ° 


jr2  r=  mx  +  c, 


we  have 


v.y, 


HT,X 


Ti  *l 


L.Xi 


L,XD=  L,y 


~S" 


v,y2 


Fig.2       Material      balance       stream: 
for     top   tray. 
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L  ,  x ,        A 


F,  xF 


~L 


HB>X; 


V,y2 


.  B  f  x. 


Fig.  3    Materiel     bolonce      streams 
for     bottom    tray. 


L(mx,  +  o)  +  V{mx„  +  e)  -  Lx.  -  V(mx,  *  o)  =.-  H,_  —i 
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dx- 


which  o>i  rearranging  yields 


dx^ 

K„,  -tr-  =  Lmx,  •;•  Lo  +  Vux2  ■;■   Vc  -  Lx,  -  Vax.  -■  Vc 


and  finally 


5X1   (Lra  -  L  -  Vm)  _  .  Vta  *  *  M  m 

__  „      __       Xi  +  __  x2  +  __  (1) 

For  the  bottom  tray,  the  light  component  balance  gives 


Input  -  Output  =  Accumulation 
I,x1  +  Fxp,  -  Vy2  -  Bx2  .  HB  ^2 


Hence 


fo-9 

HB  IF  =  LX1  *  FXF  "  V(mX2  +  0)  "  BX2 


A  total  balance  around  the  bottom  tray  shows  that 

B  «  F  +  L  -  V  (la) 

whereby 

EBinf  3  Lxl  t  FXF  "  VmX2  "  V°  "  FX2  "  LX2  *  VX2 
which  on  rearranging  yields 

dx2   Lxl    (Vi  t  F  t  L  "  V)  ,  j  PxF   Vc  ,9l 

dt  -  -fiB  "  -     Hfi    —  X2  +  1B  '  H~ 

The  overall  flow  rate  balance  around  the  entire  system  gives 
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P  =  D  +  B  ( 3 ) 

while  the  overall  light  component  balance  around  the  entire 
system,  gives 

Px„  =  DxT,  f  Bx,  (3a) 

Equations  (1)  and  (2)  describe  the  unsteady  state  behavior  Of 
the  system. 

The  following  specific  values  are  assigned  to  the  various 
parameters  of  the  system. 

F  a   0.5  lb  mole/mln. "      x  =  0,6$ 

P 

L  b  1  lb  mole/min.         H  =  1  lb  mole  (4) 

V  ■>  1.33  lb  moles/min.     H  m   2.5  lb  moles 

The  linear  vapor-liquid  equilibrium  employed  is 

y^   ,,  0.W  Xj^  +  O.56,       1  .  1,  2  (4a) 

Equations  (1)  and  (2)  have  been  used  for  simulation  on   an  EAI 
TR48  analog  computer  to  obtain  several  transient  responses  and 
phase  plane  plots.   The  analog  computer  f3ow  charts  and  scaled 
equations  are  given  in  Appendix  A.   Also  for  the  purposes  Of 
determining  the  physically  realizable  bounds  of  the  system,  the 
steady  state  and  limiting  case  analyses  of  the  system  have  been 
carried  out.   The  results  of  the  analyses  are  given  in  detaij  in 
Appendix  B. 

2,   NATURE  OF  THE  PROBLEM. 
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Using  the  parameter  values  established  in  section  1,  the 
set  of  simultaneous  equations,  equations  (1)  and  (2)  with 
the  right  hand  side  set  equal  to  zero  can  be  solved  to  obtain 
the  steady  state  values  of  x  and  x   ,  which  represent  the  liquid! 
composition  of  the  top  tray  and  the  liquid  composition  of  the 
bottom  tray  and  also  the  product  stream,  respectively.   The 
results  are 


X  m    0.6300  (5) 


and 


X      m    0.276?  .  (6) 

z 

Since 


xD  „  yx  ,-  mrx  +  c. 


Vic  have 

x    =  O.kh  x   0.6300  t  0.56  (?) 

B   0.837 

If,  at  an  instant,  the  feed  composition  x  is  instantaneously 

F 

changed  from  its  steady  state  value  of  0.65,  this  change  will 

constitute  a  disturbance  to  the  overall  system  which  in  turn 

will  eventually  result  in  a  new  steady  state. 

Let  the  initial  steady  state  (x.  =  O.63OO,  x  =  O.2767) 

corresponding  to  x     =  O.65  be  designated  by  S. ,  and  at  some 
t  -*■ 

instant  let  x  be  instantaneously  increased  to  x  =■  0.75, 
F  F 

(i.e.  we're  assuming  a  known  disturbance  to  have  occurred).   As 
a  result  of  this  disturbance,  x  will  tend  to  approach  a  new 
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steady  state  value  as  can  be  seen  in  Figure  k   by  the  path  OP. 
It  is  possible  to  determine  a  value  of  reflux  rate  L  which  will 
steer  x  from  a  certain  point  on  the  path,  say  the  point  P,  back 
to  its  initial  steady  state  value  along  FQ  asymptotically. 
However,  a  better  policy  would  be  to  change  the  reflux  L  not  at 
y   but  at  0,  and  thus  steer  x     to  its  desired  steady  state  value 
along  OKH  thereby  resulting  In  a  time  saving  of  OQ  less  OB. 
Such  a  value  of  the  reflux  rate  can  be  obtained  from  equations 
(1),  (2)  and  (3). 

From  equation  (3)  we  have 


Fx  -   DxD  +  Bx2 


=  (V  -  L)y1  +  (F  +  1  -  V)z 


=  (V  -  L)(mx,  +  c)  f  (F  +  L  -•  V)x 
3  ?. 


This  loads  to 


Fx 


F  -  (V  -  L)(mxj  +  c) 

*2  = Fim (8) 

From  equation  (1),  at  steady  state,  we  get 

V(mx,  -  mxo) 
L  ■  m'xrTc~-"xr  (9) 

and  from  equation  (2)  at  steady  state  we  obtain 


Vc  -  FXp  -  Vx,,  ■»•  VmXg  +  Fx? 


x  -  x 
1    2 


where  equations  (8)  and  (10)  contain  the  new  value  of  x  . 

F 
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By  a   trial  and  error  method,  an  L  is  first  assumed  and  x 

is  calculated  from  equation  (8).  Next  this  value  of  x     is 

2 

used  in  equations  (9)  and  (10)  and  L  is  calculated  from  each  of 

equations  (9)  and  (10)  and  checked  with  the  original  assumed 

value  of  L.  The  L  that  satisfies  equations  (8),  (9),  and  (10) 

shall  be  designated  L'.   The  value  of  L"  for  the  reference 

system  has  been  found  to  be  0.9152  lb  aoles/Mn  and  the  steady 

state  resulting  from  it  is  S  =■   (x.   =   O.6300,  x?  =_•  O.3067) 

The  increased  value  of  x_  is  consistent  with  the  material 

2 

balance  since  the  feed  composition  has  been  increased  from  0.65 

to  0.75.  It  can  be  seen  from  Figure  h   that  the  overhead 

composition  x  can  be  returned  to  it's  Initial  undisturbed 

state  in  at  least  two  different  ways,  one  of  which  (path  OKK) 

achieves  the  objective  in  a  shorter  time  than  the  other. 

Kosjever,  even  the  use  of  path  ONR  takes  a  considerable  amount  of 

tljae,  almost  to  the  extent  that  it  becomes  Impractical  to  Walt 

for  x  to  return  to  its  undisturbed  state.   Hence  a  control 
D 

policy  has  to  be  devised  for  the  manipulation  of  L  whereby  x 
is  restored  to  its  undisturbed  state  in  the  shortest  possible 
time.  This  is  what  constitutes  time  optimality. 

3.   APPLICATION  OP  THE  MAXIMUM  PRINCIPLE  ALGORITHM. 

The  Maximum  Principle  will  now  be  applied  to  the  reference 

system  so  as  to  investigate  how  the  reflux  rate  should  be 

manipulated  in  order  to  attain  the  state  S  in  the  shortest 

2 

possible  time.   A  complete  derivation  of  the  Maximum  Principle 
Algorithm  is  discussed  in  the  Appendix  C.   Stating  the  problem 
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more  explicitly  we  havei 

Given  a  system  with  the  performance  equations 

^  .  JLin  -  Vm  -_Lj      Vm      Lc  ,   . 

dt  -      HT       Xl  +  HT  X2  +  HT  U1> 

dxo   Lx-,    /„„  ,  v.  .  r  _  »\  Fx_  -  Vc 


"2 


2:    (Vm  +  F  +  L  -  V)  _   . 

,1  "  "  A^  T 


H 


dt  -  Eg         HB  2      HB 

we  wish  to  determine  L(t)  which  moves  the  system  from  the  initial 

state  at  S1  m   (x.  =  0.6300,  x  n  0.276?)  to  the  final  state  at 

S2  =  (x.  =5  O.63OO,  x?   =:  O.3067)  in   a  minimum  period  of  time.   In 

other  words  we  wish  to  determine  L(t)  so  as  to  minimize  the 

objective  function  S  where 

T 
S  «  /   dt  (12) 

o 

where  T  is  unspecified.   Physical  realizability  of  the  system 
requires  that  the  range  of  the  control  L(t)  must  be  finite 
(cf.  Appendix  B),  that  is, 

L  .   <  L  <  L  (13) 

min  •-   —  max  " 

Operation  at  the  lower  Dimit  corresponds  to  the  minimum  reflux 

L  =  L   ,  and  operation  at  the  upper  limit  corresponds  to  the 

maximum  reflux  L  =  L   . 
max 

If  we  introduce  an  additional  state  variable  such  that 

T 

x  (t)  =  ;    dt 

->       o 
we  have 


1? 

-^2  =  l,  i  (o)  =  o  (a*) 

The  Kamiltonian  function  Is  defined  as  [see  equation  (C-8)] 

it,,        .,,,   ri        Viax.,      Loz-i    T 
H  _  (Lm  -  Vm  -  L)       2     1   Jj.  _ 

H        "11  +  H   "1  +  H   *  H„  T.z2 
rp  j  qi      3 

iVm  ±F_±_L  -  Vi        Fxp  -  Vc 
B  B 


The  adjoint  differential  system  is  [see  equation  (C-9)] 

"dt  =        H        ^   fT  2  l   ' 

T  B 

dz„     Vmz,    /,,    ,.   T   ,.\ 

T  B         c 

"df=°  '  tis) 

The  boundary  conditions  on  the  adjoint  variables  are  [see 
equation  (C-9)] 

z-,  (0)   unspecified      z   (T)   unspecified 
x  1 

z2  (0)   unspecified      7.?   (T)   unspecified       (19) 
z   (0)   unspecified      e  (T)   =   1 


Equation  (18)  along  with  the  final  condition  on  z  (t)  implies 

z^  =  1  ,      0  <  t  <  T  (20) 

The  fact  that  the  final  time  T  is  unfixed  implies,  at  every 
moment 
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bz,Xt        x-,  7n        oz-,        z-i's-o       x0Zq  Vmx0z, 

ibxh   a  _    »       g  H         +  +   -g         -  /i,  + 

T  T  T  B  B  T 

Vmx1z1       Vmx2z2        ^P  "  V'X2Z2        ^FxP  "  Vc^22 
_  ._„__.  __  TT~       _.  +  _     __ 

1'        B  B  B 

+  X  =  0  (21) 

[see  equation  (A-39)]»  Also  the  Hamiltonian  has  been  written 

so  as  to  show  it  linear  in  the  control  variable  L.   Inspection 

of  the  Hamiltonian  allows  us  to  determine  the  basic  structure 

of  the  time  optimal  control  policy  as  being  a  bang-bang  policy. 

That  is,  the  control  variable  L  assumes  either  its  maximum  value 

L    or  its  minimum  value  L    as  the  system  is  transferred  from 
max  min 

an  initial  state  to  the  specified  final  state.  Of  course,  once 
the  specified  final  state  is  reached,  the  control  variable 
must  be  switched  once  again  to  L'  to  maintain  the  new  steady 
state.   The  conditions  for  the  Hamiltonian  to  be  a  minimum  are 


max      HT   "  HT  +  HT  +  Hfi     HB  ; 

„  ,mzlxl   xlzl    czl   xlz2   X2Z2,  ^  n 

Vn  "  (~h;     Hj"  +  H^  +  ~  "  "H"^  *  ° 


(21a) 


In  the  case  where  the  coefficient  of  L  (sometimes  known  as  the 
switching  function)  vanishes  over  a  finite  length  of  time  we 
have  the  possibility  of  singular  control.  In  the  likelihood  of 
singular  control  the  control  variable  may  take  on  values  which 

are  intermediate  to  I,    and  L  M    ;    hence  the  name  intermediate 

max      min 

control  is  also  used  in  place  of  singular  control. 


' 
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The  Maximum  Principle  now  requires  that  the  system 
equations,  equations  (11)  and  (14),  be  integrated  simultaneously 
with  the  adjoint  equations,  equations  (16)  and  (17)  such  that 
the  two  point  boundary  conditions 

x  (0)  ^  0.6300        x , (T)  =  0.6300 

1  1 

X  (0)  =  0.2767  X  (X)  =:  0.3067 

2  2 

x  (0)  =0.  x  (T)  unfixed  (21b) 

3  i 

Z   (0)  unfixed        z  (T)  unfixed 

1  •*- 

z  (0)  unfixed        8«(T)  unfixed 

2  2 

are  satisfied  and  assures  us  that  the  tiiae  optimal  control  is 
one  which  minimizes  the  Hamiltonian  at  every  point  of  its  response. 
Since  the  final  time  is  unspecified  this  mlnltauiu  of  the 
Hanlltonlan  is  zero,  as  indicated  in  equation  21  [also  see  equation 
CC-39)]. 

4.   COMPUTATIONAL  PROCEDURES 

The  numerical  solution  of  this  problem  Involves  guessing 
the  initial  values  of  z~    and  z  such  that  the  final  conditions 
x  (T)  and  x  (T)  may  be  reached  simultaneously  when  the  sec  of 
five  differential  equations  Is  integrated  in  a  forward  manner. 
Actually  the  guesswork  csn  be  reduced  to  the  guessing  of  only 
one  of  the  initial  values  -  either  z  (0)  or  z„(0)  -  by  makir.g 
use  of  equation  (21),   If  all  the  parameter  values  specified 
in  equations  (4)  and  (4a)  (except  L)  are  substituted  into 
equation  (21)  along  with  the  initial  values  of  x.  and  x  ,  we 
get 
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Co.2?73za(o)  -  0,63002l(0)  +  0.56z  (0)  +  0.2521z2(0) 

-  0.1108z  (0)]L(0)  +  0.16258^(0)  -  0.3696z-(O) 

-  0.065z  (0)  4   0.O923Z  (0)  -  0.l486z2(0)  +1=0 

which  on  rearranging  yields 

i 

(0.20?1L(0)  ■-  0.2071)2^0)  -i-  (0.Hfl3L(0) 

-  0.1213)z  (0)  +  1  a  0 


1  +  [0.141 31,(0)  -  0.1213]z2(0) 

z1(o)  m  -      67WnnW~-~o'72ofi'~'  *22* 


This  equation  indicates  that  it  suffices  to  guess  the  initial 

value 

1.(0). 


value  of  z„  only,  since  z  (0)  novr  is  dependent  on  z  (0)  and 


The  system  equations,  equations  (11),  (16)  and  (1?)  have 
been  Integrated  on  an  IBM  system  360  computer  with  the  help  of 
the  IBM  supplied  subroutine  HKGS  (Runge-Kutta  Gill  subroutine). 
The  HKGS  subroutine  solves  a  set  of  simultaneous  differential 
equations  with  given  initial  conditions.  The  overall  logic  flow 
diagram  for  the  trial  and  error  solution  can  be  seen  in  Figure 
5  where  the  following  procedure  is  indicated  in  the  corresponding 
numbered  boxes. 

(1)  The  known  initial  values  namely,  x, (0),  x  (0)  and  x  (0) 

12        3 

are  read  In. 

(2)  A  choice  of  either  L    or  L    is  made  for  L. 

max     min 

(3)  The  initial  .value  of  z     Is  assumed. 

{k)     The  corresponding  initial  value  of  2  is  calculated  from 
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Sel 


L=L 


.:,:;■; 


Assume 

2,(0) 


Calculate 

z^Olfrom 

equation  (22) 


9- 


Calculate  value 

of 

switching 

function 


.-•■.-- 


10 


It.   si?:-,  of  switch- 
ing function  cons- 
istent with  choice 
of     L? 


No 


Select  L 
according  to  sijr 
of   switching 

function   in  boxll 


Ye  ■ 


1 

READ  ~]         s~ N 

0),x,(0)/— -  (      START       ) 


•:. 


Sat  total 
'  integrating 
time  at  V 


Call 
RK6S 


Set 
At  =  0-OI 


-£. 


jt  =  t  +0  01 


Fig. 5.   Computer    logic   diagram  for  trial  and  error 
solution    of  equations  (II),(I4),{!6)  and  (17). 
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equation  (22). 

(5)  1'he  integration  procedure  is  set  to  terminate  at  some  large 
value  T' . 

(6)  The  increment  in  time  is  set  at  0.01. 

(?)  Subroutine  HKGS  is  called,  and  the  integration  commences. 
HeD'e  mention  must  be  made  that  the  commands  to  be  executed 
in  boxes  8  to  \Ur   inclusive,  take  place  over  each  t. 

(8)  The  current  values  of  x  ,  x  ,  x  ,  z  t  z  and  t  are  written 

out . 

(9)  The  value  of  the  switching  function  is  calculated, 

(10)  The  sign  of  the  switching  function  is  checked  for  consistency 
with  the  choice  of  L  in  box  2.   If  the  choice  of  L  is  con- 
sistent with  current  sign,  then  command  transfers  to  box  12, 
if  not,  to  box  11. 

(11)  A  choice  of  L  is  made  to  match  the  current  sign  of  the 
switching  function. 

(12)  The  current  value  of  the  Hamiltonian  is  written  out.   The 
reason  for  doing  this  is  to  check  whethei*  H  maintains  a 
constant  value  of  zero  during  the  transient  response. 

(13)  The  time  is  incremented  to  the  next  higher  value. 

(,1k)   The  current  value  of  time  is  compared  with  the  total 
Integrating  time  specified  in  box  5> 
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For  the  initial  guess  zJ.O)   to  be  the  right  one,  the  desired 

final  values  of  x.  and  s  Should  oocur  at  the  same/lt.  Also  it 

is  Imperative  that  the  Bamlltonian  vanish  throughout  the  optimal 

transient  response. 

The  above  procedure  of  calculation  was  applied  to  the 

reference  system  after  selecting  suitable  values  for  L    and 

max 

L  ,  .   Initially  a  preliminary  case  was  Investigated  wherein 
nin 

L    was  set  at  0.930  and  L    at  O.896.  The  optimal  policy  was 
max  mi  n 

found  to  be  to  start  with  L    and  then  to  switch  to  L  ,   after 

max  lain 

2. 46  minutes.   The  final  state  was  attained  via  L  ,   in  an 

nun 

additional  1.28  minutes.   For  the  preliminary  case  above,  the 

values  of  L    and  L  ,   were  selected  in  an  arbitrary  manner; 
max     min 

however,  for  furthur  analysis  of  the  system,  L    and  L  ,   were 

max      mln 

selected  in  a  manner  symmetric  about  L'  =  0.9152,  (which  is  the 

reflux  rate  that  maintains  the  desired  final  state)  as  seen  in 

Figure  6.   A  point  symmetric  with  L  =  0,833  on  the  right  side 

of  L'  was  found,  to  be  L  ^   0.997'u  Hence,  L  =  0,833  and 

L  ■  0.997*}  were  taken  as  the  0%   and  100$  points  of  the  range 

within  which  L    and  L  .   were  selected.   (The  maximum  and 
max     min 

minimum  bounds  en  L  have  been  established  in  Appendix  B).   The 
different  cases  investigated  weret 

Case  1  L  ,   =  O.8330  L         m   0.9974 
mm  max 

Case  2  L  ,   =  0. 86588  L    -   0.96452 
min  max 

Case  3  L  ,  =  0. 89876  L    =  0.93164 
min  max 


In  each  of  the  above  cases  the  average  of  L    and  L    is  L1. 

max      min 


'A 


0.833  0.&I52  Q9974  |333 

Case  I        Lmax=  0.99740  100% 

Lmin  =0.83300  0% 

Ca<,e2       Lmax  =  0.96452  80% 

Lmin  =0.86588  20% 

Co,e3       Lmax  =0.93164  60% 

Lmin  =0.89876  40% 

Case  4      Lmax  =  1-3330     ' 
Lmin  "0.8330 


Fig.  6.    Physical    bounds    within    which    L 

is   constrained    for    cases  I,  2, 3  and 4. 
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In  addition  to  the  above  tfrree  cares  a  fourth  case  was  Investigated 
wherein  the  whole  horizontal  line  in  Figure  6   was  included.  For 
this  case  we  have 

Case*   LBin,-  0.833   ^  -  1.333 
Apart  from  the  preliminary  case  where  L    «  0,980  and 

TdcX'K. 

L  ,     =  0,896,  Case  1  was  also  solved  by  the  trial  and  error 
nun 

method  outlined  in  Figure  5,  wherein  the  adjoint  equations  were 
solved  along  with  the  system  equations.  However,  this  method 
was  found  to  be  tedious  and  time  consuming  and  in  place  of  It  the 
phase  plane  method,  used  by  Siebenthal  and  Arls  [8],  was 
employed  to  obtain  the  solutions  for  Cases  2,  3:  and  k.      In  using 
the  phase  plane  method  Pontryagin's  Maximum  Principle  was 
employed  only  to  find  that  bang-bang  control  was  necessary  and 
then  the  adjoint  equations  were  discarded  and  phase  plane 
diagrams  were  used  to  obtain  the  solutions.   This  enables  us  to 
circumvent  to  some  extent,  the  inherent  difficulties  of  the  two 
point  boundary  value  problem.  Before  going  into  the  description 
of  the  phase  plane  method,  the  following  two  facts  are  worth 
noting. 

1.  The  state  equations  are  autonomous,  (i.e.  time  does  not 
appear  explicitly  on  the  right  hand  side  of  the  equations) 
Theorems  on  the  uniqueness  of  solutions  for  autonomous  systems 
[10]  guarantee  that  for  a  given  value  of  the  control  variable 
L,  trajectories  in  the  phase  plane  do  not  intersect  at  any  point 
except  at  the  steady  state  point  corresponding  to  L. 
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2-     It  Is  known  that  during  the  transient  period,  the  optimal 

control  policy  is  bang-bang-  and  furthermore,  the  number  of 

switches  In  the  value  of  the  control  variable  will  be  kept  at  a 

minimum  since  the  time  taken  to  get  from  the  initial  state  to 

the  final  state  would  increase  as  the  number  of  switches  is 

increased  [ll]. 

In  Figures  7,  8,  and  9  phase  planes  are  presented  for  the 

preliminary  case  where  the  values  of  L  .   =  0.896,  L'  =  0.9152 

min 

and  L    a  O.98O  respectively,  (this  is  the  preliminary  case 
max 

already  solved  by  the  trial  and  error  method  on  the  digital 

computer).  .In  Figure  7,  point  A  is  the  steady  state  point,  and 

in  Figure  9,  point  B  is  the  steady  state  point.   The  steady 

state  point  S  on  Figure  8  is  the  desired  final  state.   This 

point  is  also  shown  on  Figures  7  and  9.      Since  the  system  is 

autonomous,  only  one  curve  passes  through  the  point  S  In  the 

Figures  7  and  9  [10]  these  curves  are  lettered  DS  and  E3 

respectively.   If  Figures  7  and  9  are  superimposed  as  shown  in 

Figure  10,  it  becomes  obvious  that  there  are  an  infinite  number 

of  paths  consisting  of  alternating  L  ,   and  L    response 

nun     max 

segments  which  connect  a  given  initial  point  to  the  point  S  . 

However,  all  of  these  possible  paths  approach  S  via  either  the 

I.    segment  E5„  or  the  I.  .   segment  D3  . 
max  2        min  2 

Closer  inspection  of  Figure  10  shows  that  it  is  possible 

to  connect  every  point  on  the  x  ,  x?  plane  to  the  point  S  by 

means  of  a  path  consisting  of  one  I.    and  one  L    segment  or 

max         min 

vice  versa.   Such  a  path  is  S  FS  in  Figure  10  where  S 
corresponds  to  the  initial  state  of  the  reference  system 
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Fig,  7. 
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Phase    plane    of    reference    system 

with    L  =  0-896  ib  moles/min.  (Preliminary  Case). 
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Fig.  8.    Phase    Plane    of    Reference    System 
with  l>  0.9152      lb.  moles /mia. 
Preliminary  Cose. 
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Fig.  9.    Phase   plane    of  reference   system 
with    Lmax  =  0.980    lb.mo!es/min, 
(Preliminary  Case). 
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Fig.  10 .    Superimposifion        of      Figures    7    end  9 
{Preliminary    Ccse)- 
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(0.6300,  0.2?67).  This  path  has  only  one  switch  at  point  F, 

Also  the  path  S  FS  oonf litis  the  phase  plane  method  of  solution 

since  the  point  P  agrees  with  the  point  in  the  digital  computer 

solution,  previously  obtained,  at  Which  the  snitching  function 

changes  sign  and  the  control  policy  switches  from  L    to  L  .  . 

max     mj.n . 

Hence,  the  only  control  policies  which  satisfy  the  Maximum 
Principle  equations  contain  just  one  control  switch  which  occurs 
when  the  response  Intersects  the  curve  DS  E  in  Figure  10.  Curve 
DS  E  is  called  the  "switching  boundary"  due  to  the  fact  that  it 

divides  the  phase  plane  into  regions  of  L    and  I.    operation. 

*  *  max      min 

In  Figure  10,  operation  in  the  region  above  DS^E  is  with  L 

while  operation  in  the  region  below  DSgE  is  with  Lnax>   This 

method  of  phase  plane  analysis  was  also  used  in  obtaining  the 

optimal  control  policy  for  Cases  1,  2,  3>  and  &  as  follows! 

First  the  phase  plane  plots  were  obtained  on  an  analog 

computer  by  integrating  the  system  equations  only  and  then 

phase  plane  plots  similar  to  Figure  10  were  obtained  individually 

for  the  4  cases  by  superimposition  of  L  .   and  L  _  response 

mln      max 

phase  planes.  These  can  be  seen  in  Figures  11,  12,  13  and  Ik. 

Next,  on  the  digital  computer,  the  system  equations  were  first 

integrated  forward  in  time  from  the  point  S,  using  L    as  the 

control  variable,  and  next  integrated  backward  in  time  from  the 

point  S  ,  using  L    as  the  control  variable.   By  comparing  the 
2        mln 

computer  outputs  the  common  point  F  was  located.   Once  more  a 

forward  integration  was  carried  out  from  S  and  using  L    as 

1  -max 

control  variable.  Also  this  time,  a  command  was  Included  for 

the  control  variable  to  switch  from  L„ „„  to  L  .   when  F  was 

max     mln 
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Fig. II.  Superimposi.fion  of  LmcixQria  t-min  response  phase  planes  for 
Case  I.  Paths  directed  to  point  B  are  l.-^*  response  curves 
orxi    paths    directed    to    point  A  ere   Llnin   response     curves. 
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Fig.  12.      Supe (imposition     of    Lmax    and  Lmjn  response    phase    planes 

for    Case-2.   Paths  directed    to  point    B    ore    Lmox    response 

curves    ond    paths  directed    1o  point    A    are    Lmjn     response 
curves. 
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Rg.13.  SuperimposMon     of  LmQX    and    |_mjn     response  phase    planes 

for   Case  3.    Paths  directed     to    poirrt   B    ore  Lmox  response 

curves    ond     paths  directed     to    point   A  ere  Lft,jn    resp 
curves. 
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Fig  14.   SuperimposWon  of    Lmm    end    L,njn    response    phase   planes  for 

Case  4.     Paths    directed     to    point    B  ore    LmQX    response    curves 
end     paths    directed     to    poinl    A    ore    Lrnjn    response    curves. 
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reached.  The  integration  was  then  allowed  to  proceed  with 

L    until  S  was  reached  at  which  stage  the  computation  vias 
Kiln       <- 

terminated.  The  value  of  x  at  the  final  state  gave  the  desired 

minimum  of  the  objective  function. 

Table  1  shows  the  results  of  the  above  computations  as 

applied  to  Cases  1,  2,  3  and  4.  In  each  case  control  starts 

with  L    and  switches  to  L  .   at  the  switching  time  t  . 
max  mln  s 

Figures  15,  16,  1?,  and  18  show  the  transient  response  of  the 
system  and  the  corresponding  control  policy  for  Cases  1,  2,  3  and 
4. 
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TABLE  1 


RESULTS  OF  REFERENCE  SYSTEM  ANALYSIS  FOR  THE  TIME 
OPTIMA],  PROBLEM 


Order  in  which       Switching  time    Final  time 

T,  inin 


Case  1     L    ■  0.9974  2.24  2.80 


Case  2     L     0.964-52  3.41  4-. 05 


control 

is  applied 

t  ,  mln. 
s 

L 

max 

0.9974 

2.24 

L  ,     = 
mln 

0.833 

L    = 

max 

0.96452 

3.41 

L  ,   = 

mln 

O.86588 

L 

max 

0.93164 

'   7.63 

L  ,  = 
min 

0.89876 

L       B 

max 

1.333 

0.56 

L  ,  = 
min 

0.833 

- 

Case  3     L     0.93164      '   7.63  8. 16 


Case  4     L    =  I.333  0.56  1.68 


0.64 
< 

y^ 

-"\ 

0.63 

, ' 

\ 

0.33 

£J    0.30 

- 

—■ — "    - 

0.H7 

0.9974 

J      0.3i52 

_ 

0.833 

- 

-- 

Time  miniutos 


Fig- 15.  Responses    of  reference    system  to  bang- 
bong    policy  (cose  I)  for  time    optimal   pro- 
blem   (path  S,FS>of  Fig. 1 1.) 
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Fig.  16.  Response  of  reference  system  to  bang- 
bang  po!icy(Case.2.)  for  time  optimal  pr 
oblem(path  S.FS.of  Fig.  12.) 
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F"ig.'i7      Response     of     reference     system    to    bang- 
bong     policy  (Case  3)    for    time     optimal 
problem  (path  S,FS2  of    Fig. 13). 
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Fig.  18.     Response      of     reference     system      to     bong-bang      policy 
(Case. 4.)     for     time      optimal      problem.  (Path  S.FS^of 
Fig. 14.). 


2:0 


42 


CHAPTER  3 


SYSTEM  WITH  A  TOP  TRAY  DESCRIBED  BY  THE  MIXING  POOL 
MODEL  WITH  TWO  TANKS  IN  SERIES 


In  this  chapter  we  shall  consider  the  same  problem  as  in 
the  previous  chapter,  i.e.,  the  time  optimal  problem,  except 
that  the  top  tray  of  the  column  will  be  described  by  the  mixing 
pool  model  with  two  tanks  in  series j  henceforth  referred  to  as 
the  TTIS  system. 

1.   DEVELOPMENT  OF  PERFORMANCE  EQUATIONS 

The  following  assumptions  and  column  characteristics  will 
be  appended  to  those  already  made  in  Chapter  2. 

(a)  The  top  tray  is  described  by  a  mixing  pool  model  with  two 
completely  mixed  pools  of  equal  volume  in  series. 

(b)  The  vapor  rising  from  the  bottom  tray  is  assumed  to  divide 
in  half  before  entering  each  pool  of  the  top  tray  and  likewise 
after  emerging  from  each  pool,  the  vapor  streams  merge  into  a 
single  stream  before  being  totally  condensed. 

(c)  The  overhead  liquid  composition  is  equal  to  the  mean  of  the 
compositions  of  the  vapor  streams  emerging  from  the  two  pools 
I.e. 

XD=l(yi+y2)  <23) 

Figure  19  shows  a  block  diagram  of  the  mixing  pool  model  system 
with  the  various  flow  streams  involved. 

The  dynamic  behavior  of  the  column  in  the  transient  state 


^3 


Fig.  19.     Distillation     column     with    top    tray    described  by 
mixing    pool     model     with     two    tanks     in   series. 


H 


may  now  be  written  in  the  different  rial  form  by  performing  a 
mass  balance  on  each  trr.y  along  with  overall  material  and  light 
component  balances,  as  follows  .(see  Figures  20,  21,  and  22); 
For  the  top  tray  pool  1, 

Input  -  Output  =  Accumulation 

Lx  +  Jf  y  -  Lx  -  V  y  mh^k 

D   2  *3     1   2  1    2   at 


Substituting   from  equation   (23)    gives 

K     dy. 

|(y2  +  y2)  +  \  jf3  -  lXi  -  |  yx  =  -J  -^ 


Also  since 


y     =  mx     +  c    ,  i  =  1,    2, 


we  have 


|(mx1  +   c  +  mx£  +   c)  +  f(mx3  +   c)   -  Lxj  -  jdBXj  +   c) 


2    at 


and  finally  rearranging  gives 


dx- 

ir  -  ■      -e~-      -1 T  ht-2 


1  _  .(I.m„ -  2L  - Vml      Lm     Vm     2Lc         ,  ,  . 

_  _  J_ __ :J-  x..  +  Tf-X„  +  „-X„  +   „  IZ*J 


For  the  top  tray  pool   2, 


Lx1  +  fy3  -  Lx2  -p-2  =  T7[r 


Since 


L.x. 
t •  i 


^ 


V,y, 

2    v 


L,x, 


^5 


Fig.  20.     Material     balance    streams    for 
poo!  S      of    the     top    tray. 
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Fig. 21      Material     balance    streams    for 
poo!  2     of    the    top    tray. 
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Fig. 22      Material     balance    streams    for    the 
bottom    tray. 


yx  =  mxj  +  c  ,      1  =  1,  2,  3. 


wc  have 


^  *  i<-3  ♦•>-*,  -*-*♦•>-££ 


which  on  rearrangment  yields 


And  for  the  bottom  tray, 

Fxp+  Lx2  -  Vy3  -  Bx^  HB  -^ 


Since 


y  =  mx„  +  c, 


we,  have 


k8 


dx2   2L  .,    (Vm  +  2L)      Vin  ,„r, 

"dF  =  BT  Xl  "  J— iL™   x2  +  H~x3  (25) 


dx_ 

Fx     +  Lx„   -   V(mx„  +   c)    -  Bx_  =   Ht,  —-2 
F  2  J  3  B    it 


and  reari'anging  yields 

^2  ,  i  x    -  fo»  +  y  +  l  -  v).    .  2e  .  Ic  .  6 

dt    " '   H        2  H  x3  +    H  H  UbJ 

B  B  B  B 

Equations  {2k),    (25),  and  (26)  are  the  performance  equations  of 
the  TTIS  system.  The  various  parameter  values  are  identical  to 
those  in  equation  CO  and  once  more  the  linear  vapor-liquid 
equilibrium  i.e. 


y  =  O.Wtx    +  0,56  ,      i  =  l.  2,  3 


1*9 


Is  obeyed.  The  above  performance  equations  were  used  for 
simulation  purpose  on  the  analog  computer  whereby  the  various 
phase  plane  plots  and  transient  responses  were  obtained.  The 
analog  computer  circuit  diagrams  are  explained  in  Appendix  A. 
Also  the  steady  state  and  limiting  case  analyses  are  performed 
in  Appendix  B. 

2.   NATURE  OF  THE  PROBLEM 

As  in  the  case  of  the  reference  system,  equations  (2k),    (25), 

and  (26)  were  solved  simultaneously  (after  setting  the  left  hand 

sides  equal  to  zero)  to  obtained  the  steady  state  values  of  x  , 

x  and  x„.   These  were  found  to  be 
2      3 

x  =  0.7146? 


and 


x  =r  0.60909 
2 


x  =  0.24899 


Since 

x  =  ^(mx  +  mx  +  2c) 

=  0.22  x.  +  0.22  x,  +  c,  (27) 

we  have 

x  =  0.22  x  0.7146?  +  0.22  x  0.60909  *  O.56 
=  0.85123. 


Let  the  Initial  steady  state  corresponding  to  x  =0.65 

F 

be  designated  S«  2  (xn  =  .71467,  x„  =  .24899).   If  x  is 
11  2  p 

instantaneously  increased  to  x  =  0,75  the  transient  response 

F 

given  by  curve  OP  in  Figure  23  results,  showing  the  displacement 
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of  x  at  S'  to  a  higher  steady  state.  Once  more  a  value  L'  of 

the  reflux,  rate  can  be  corcDuted,  which  villi  return  x  to  its 

D 

initial  steady  state  either  via  PQ  or  ONE.   Such  an  L'  can  be 
deternlned  as  follwsj 

From  the  overall  light  component  balance,  we  have 

FxF  a  taD  +  Bx3 


Fxp  -  DxD 

x3  "    35  ~ 


FxF  -  (V-L)xD 

P  +  L  -  V 


(28) 


At  steady  state  equation  (24)  gives  rise  to, 

Vc  ~  Fxv,  +  Vmx-j  -  Vx-j  +  Fx-, 

L  =  1— 2 2 1.  (29) 

x„  -  x„ 

2    3 
equation  (25)  to 


Vmx?  -  Vmx. 
2(x1  -  x~T 


L  =  __2_^^  ,  (30) 


and  equation  (26)  to 


Vmx-,  -  Vbjx-5 


x,  -  BjS  +  mxp  +  2c 


(3D 


The  trial  and  error  procedure  for  deterralng  L  from  these  four 
equation.^  is  similar  to  the  one  used  in  the  case  of  the  reference 
system  i.e.  a  value  for  L  is  assumed  e.nd  x  is  calculated  from 
equation  (20).  Next  this  value  of  x  is  used  and  L  is  calculated 
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from  each  of  equations  (29).  (30),  and  (31)  and  checked  with 

the  assumed  value  of  L.   That  value  of  L  which  satisfies 

equations  (28),  (29),  (30),  and  (3D  simultaneously  Is  designated 

L'.  The  value  of  L*  for  the  TTIS  system  was  found  to  be  0.922 

lb.  moles/min.,  and  the  steady  state  resulting  from  L'  is  S*  s 

(x  =  0.71399,  x,  =  0.60986,  x  =  0.2824?).   Even  though  L'  does 

return  x  to  its  initial  steady  state,  the  time  taken  to  do  so  is 

impractical  and  hence  an  optimal  policy  for  L  has  to  be  determined 

so  as  to  restore  x„  in  a  time  optimal  manner.  Hotel  x_  as 
D  D 

defined  by  equation  (2?)  Is  the  same  in  S'  and  S'  even  though  x 

.1       c.  X 

and  x  in  the  two  states  are'  different. 
2 

3.   APPLICATION  OF  THE  MAXIMUM  PRINCIPLE  ALGORITHM 

Once  again  the  optimum  reflux  policy  shall  be  determined  via 
the  Maximum  Principle.   The  performance  equations  of  the  system 
are  equations  (2*0,  (25),  and  {26)   and  are  repeated  below  for 
convenience. 

dXl    (Lm  -  Vm  -  2L)  _    Lin     Vm     2Lc 
dt  "       HT       xl  +  HTX2  +  H~X3  *  "K^ 

ff2_2L    -  J2L  •;•  Vm),      Vm  ,   , 

dt  -  H  xl       H     x2  +  H   x3  *3Z' 

T  T  T 

51  -  HE  H  Hfi        X3  +  ~-H-~- 

The  system  is  to  be  transferred  from  an  Initial  state  S'  5 
U1   "  °-7^(>7,   x?   =  0.60909.  x  =  0,24899)  to  the  terminal  state 
S^  B  (^  o  O.71399,  x2  =  0.60986,  x  =  0.2824?)  by  manipulating 
the  reflux  rate  L  and  simultaneously  minimizing  the  objective 
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function 


T 

s  -  /  at  (33) 

o 


where  T   is  unspecified.   The  control  variable  is  constrained  by 

L  .   <  L  <  L  (34) 

ram  -   -  max 

An  additional  state  variable  is  defined  such  that 

T 
x  (T)  =  S  a  I     dt  (35) 

4  o 


thus 


dt  "a 


According  to  equation  (C-8),  the  Hamiltonian  for  this  system 
becomes 

F  =  lLm_-  Vm-  2L).      +  La      +  |m     +  21^ 


ii 


a 


11  T  HT  A2"l  v  Hrj?A3"l  T  H  "1 


+  uj\   Z  Hm    X2Z2  +  H„X3  2  +  H 


r  yz  -   3^2*3 


(V.  t  ,+  L  -  V)     +  (ISLLl^  +  (36) 

B  B 


and  the  adjoint  differential  system,  according  to  equation 
(C-9),  becomes 

dZl  _  _  .(Lin^-Vm  -  2L)_   _  _2L 


■,  -  -r-"  *o  (37) 

T 


dt  ■  "1    HT  "2 


at 

■--. 

- 

> 

21 

+ 

i£ 

d  z 
dt 

= 

- 

Viu 

Zl 

- 

Via 

dzj, 

_ 

0 

5'+ 


I^AJSslz,  .  JLZ  (33) 

HT     2   HB  j 


+  (VtL^JLtJi^J/)  «         (39) 


(to) 


(to) 


hb      "3 


at 

The  boundary  conditions  on  the  adjoint  equations  are 


z  (0)  unspecified  z   (T)  unspecified 

z  (0)  unspecified  z.(T)  unspecified 

e  (0)  unspecified.  z  (T)  unspecified 

z.  (0)'  unspecified  zj,(T)  -  J 


From  the  final  condition  on  z^(t)  and  equation  (40),  we  obtain 
z.  »  1,      0  <  t  <  T   '  (42) 


Since  the  final  tine  T  Is  unspecified  and  the  top  tray- 
hold  up  H  ,  is  unity,  the  minimum  of  the  Hamiltonian  may  be 
rewritten  as 


min  H  =  (in-.z,  -  2x,  z,  +  mx0z,  +  2czn  +  2x,  z„  -  2x„z„ 
11      11      c   1       1      12      <c  <: 


X.jZ-.  X0Z3 

+  __^J  „  _^)L  +  vte^  .  vm^Zj  ,.  Vta2z2  -  Vaa3z2 


3E  HB 

Vmx-,2.,        /  D        ,m  Fx„  -   Vc 

_J2_2  „  <f  -  v>  +  (_JL___„,      +  x  .  0        (43j 


3  3  H  3 

B  B        J   J  B 


Once  more  the  Hamilfconian  is   shown  linear  in  the   control 


variable  L  which  implies  that  the  optimal  control  poDiey  Is 
probably  s.   bang-bang  policy  as  in  the  case  of  the  reference 
system.  The  control  policy  will  be  governed  by  the  sign  of  the 
switching  function,  i.e. 


L  -  Lmaxif  (BXlzi  "  2xlzl+  mx2Zl+  2cSl*2xlZ2-  2*^ 


HB     HB 


L  =  L   ,    if   (mx-Z--  2x.,z..-fr.  n!x0z..(.   2cz,+   2x-z„-  2x_z» 
min  11  11  ^  i-  1  12  ii  <L 


W) 


+  -|-3  -  -3-3)  >  o 


B 


B 


Now  it  is  required  that  the  system  of  equations  (32)  and 

(35)  be  integrated  simultaneously  with  the  adjoint  system, 

t 
equations  (37),  (38)  and  (39).  such  that  the  two  point  boundary 

conditions 


x  (0)  =  0.7146? 

x  (0)  ^  0.60909 

x  (0)  =  0.24899 

x^O)  =,0.0 

z  (0)  =  unspecified 

z   (0)  =  unspecified 

z_(0)  =  unspecified 


x  (T)  »  0.71399 

x  (T)  =  0.60986 
2 

x.  (T)  =  0.28247 

x  (T)  unspecified 

z  (T)  unspecified 

z  (T)  unspecified 

z  (T)  unspecified 


are  satisfied  and  the  resulting  control  minimizes  tne  Harailtonia.n 
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at  every  point  of  Its  response.  Once  again  the  final  time  is 
unspecified  and  hence  the  minimum  value  of  the  Hamlltonian  is 
zero  throughout  the  optimal  respor.ee . 

k.      COMPOTATIONAL  KHOCEDURES 

The  direct  numerical  solution  of  the  above  set  of  differential 
equations  is  considerably  more  difficult  than  that  of  the 
reference  system  since  the  present  solution  involves  the  guessing 
of  two  of  three  initial  values  z  (0),  z  (0)  and  z  (0)  instead  of 

X  c  J 

a  single  initial  value  as  In  the  case  of  equation  (22).   By 
employing  the  property  that  the  Hamlltonian  must  vanish  at  very 
point  of  the  optimal  transient  response  we  obtain  at  the  initial 
point,  an  expression  of  the  form 

.  ,     rl+(0.2112L(0)~0.2112)z-(0)+<Hl^lL(0)-0.12to)zo(0)., 

z  t  o  j  =  -  l — • fi — ' J 

0.2?32L(0)  -  0.2732 

145) 

In  order  to  eliminate  the  complexities  of  the  two  point  boundary 
value  problem,  the  adjoint  equations  were  discarded  and  once  more 
the  phase  plane  method  of  analysis  was  resorted  to.  However, 
due  to  the  three  dimensional  nature  of  the  problem  we  are  confrontec 
with  a  phase  space  rather  than  a  phase  plane.   Hence,  a  single 
phase  space  will  be  represented  by  a  set  of  three  phase  plane 
diagrams [  that  Is,  the  phase  space  of  x  x  x     will  be  spanned  by 
the  phase  planes  X-jXg,  r?x   ,   and  x_x,.   Also  the  constraints  on 
L,  i.e.  L    and.  L    are  picked  symmetric  with  respect  to  the 
reflux  rate  which  restores  the  system  to  the  desired  final  state, 


l_        . 

L 

1 i      

0-833 

0922             l-OII 

i-; 

Case  1 

LmQx=  io  no 

100  % 

Lmin  =  0-8330 

0  % 

Case  2 

L-max  =0  9754 

80% 

Lmin  h  0-8686 

20% 

Case  3 

LmQX  -  0  9398 

60% 

Lmin   =0-9042 

40% 

Case  4 

LmQX  =  1-3330 

Lrnin    =0-8330 

I'iS 


Fig.24.     Physical    bounds    within    which    L    is   constrained 
for   Cases   1,2.3  and  4  . 
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i.o.  symmetric  to  L  -  0.922  (which  was  earlier  denoted  as  L' 
from  equations  (28),  (29).  (30)  and  (3D.  This  gives  rise  to 
the  following  three  cases  (see  Figure  2>4), 


Case  X    I,    =  0.833   L    a  1.011 

o i  i  max 


Case  2    L      0.8686  L    .-_■  0.975'* 
mi  n  max 


Case  3    L  .  ■  0.90^2  h         =  0.9393 
ffl-in  max 


In  addition  to  the  above  three  cases  a  fourth  case  was 
investigated  wherein  the  entire  range  of  L,  as  determined  in 
Appendix  B,  was  considered.  This  corresponds  to  the  total  length 
of  the  horizontal  line  in  Figure  2'!.   For  this  case  we  have 

Case  h         L     0.833   L    -  1.333 
min  max 

Please  plane  diagrams  for  the  above  four  cases  were  plotted 
with  the  help  of  an  analog  computer  and  are  of  the  type  shown  In 
Figures  25,  26,  and  27.   Which  are  the  phase  planes  for  Case  k 
In  each  of  the  phase  planes  the  initial  state  is  denoted  by  S' 
and  the  desired  final  state  is  denoted  by  S'. 

Before  going  into  the  search  for  the  optimal  control  policy, 

the  above  four  cases  were  simulated  on  the  analog  computer  and 

the  resulting  phase  plane  plots  were  examined.  An  Interesting 

observation  was  made  from  studying  these  phase  planes  closely. 

Starting  from  the  initial  state  at  S" ,  the  system  travels  along 

path  S'T  or  S'Q  depending  on  whether  the  Initial  control  action 

Is  L    or  L    respectively.   Once  the  system  travels  along  one 

of  these  initial  paths,  then  no  matter  what  combination  of  L 

'  max 
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0.85 


Fig.  25.     Superimposed      phase     planes  (x,,x3)  of 
two  tanks   in  series  system    for    Case  4 . 
(Lmax.  =  1-333,  Lmin.=0-833). 
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Fig.26.  Superimposed    phase    planes  (x2,x3)  of 
two  tanks  in   series    system  for  Case  4 
(Lmin=  0-833,  Lmax  =  1-333). 
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Fig.27  Superimposed     phase    planes  (x,  ,x?)  of 
two  tanks  in  series  system  for  Case  4 
1 1^^=0-833,  LrflflV -  1-333). 
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and  L    paths  is  used,  the  desired  final  state  at  S'  can  never 
miti  2 

bs  attained.  It  is  essential  at  this  point  to  explain  the 
fact  that  even  though  it  may  appear  from  Figures  2$,    26  and  27 
that  a  path  does  exist  which  passes  through  S',  this  is  not  so. 
In  order  to  elaborate  a  little  on  this  point  let  us  consider 
Case  k  whose  phase  space  is  represented  by  the  phase  planes  in 
Figures  255  ?-6  and  27.  In  Figure  ?5r  for  instance  t  it  seems 

obvious  that  an  L  ,   path  can  pass  through  S',  however,  on  close 

nun  2 

examination  it  is  found  that  at  S'  only  the  final  values  of  x_ 

and  x^  are  satisfied  while  x  has  a  value  of  0.61W6  which  is 
3  2 

higher  than  the  desired  value  of  0.60986.  Similarly  in  Figure 

26  at  S'  the  final  values  of  x  and  x  are  satisfied  while  x  is 

found  to  be  0. 71033  which  is  lower  than  the  desired  value  of 

0.71399.  And  likewise  in  Figure  27  at  S'  the  final  values  of 

x  and  x  are  satisfied  while  x  is  found  to  be  considerably 

J      2  J 

lo.sr  than  the  desired  value  of  0.282^'7. 

The  above  piece  of  observation  leads  to  the  conclusion  that 
bang-bang  control  alone  cannot  achieve  the  desired  terminal 
state  In  a  finite  number  of  switches.   It  may  be  that  the 
terminal  state  can  be  attained  in  an  infinite,  number  of  bang- 
bang  switches  but  this  being  highly  impractical,  the  possibility 
can  be  safely  abandoned.   Therefore,  the  likelihood  of  singular 
or  intermediate  control  becomes  a  certainty  since  it  seems 
logical  that  if  the  final  state  cannot  be  arrived  at  by  purely 

bang-bang  control,  some  control  variable  other  than  L    or 

max 

L         has  to  be  used  in  conjunction  with  L    and  L   . 
Bin  max      min 

In  general,  a  singular  control  is  indicated  whenever  the 
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switching  function  becomes  identically  zero  over  one  or  more 
positive  intervals  of  t.V.uc.  When  such  a  situation  arises  the 
Hamlltonian  becomes  independent  of  the  control  variable  and 
consequently  the  Maximum  Principle  falls  to  yield  a  well  defined 
control  policy.  At  present  thei'e  is  no  generally  applicable 
analytic  method  by  which  one  can  determine,  a  priori,  whether 
a  singular  candidate  actually  represents  an  optimal  trajectory, 
except  in  the  case  of  certain  restricted  classes  of  problems 
[12].  Also  whenever  the  singular  control  is  used,  the  existence 
theorems  of  Marcus  and  Lee  [13]  does  not  guarantee  the  existence 
of  an  optimal  control  policy.  A  statement  of  these  theorems  is 
given  in  Appendix  C. 

Thus  the  problem  of  synthesizing  the  optimal  control  becomes 
one   of  determining  the  optimal  switching  hypersurfa.ee  S  which  is 
the  union  of  a  bang-bang  switching  surface  S  and  a  singular 
control  surface  S  .  As  an  optimal  trajectory  is  traced  out  in 
the  x,  t  space,  the  representative  point  penetrates  the  surface 
S  at  each  switch  of  the  bang-bang  control  given  by  equation  (44), 
When  the  optimal  control  becomes  singular,  the  representative 
point  impinges  upon  and  continues  to. follow  along  the  surface 
8  .  A  typical  situation  is  illustrated  in  Figure  28. 

In  general,  the  singular  trajectories  become  unstable 

solutions  of  the  canonical  equations,  equations  (32),  (35),  (37), 

(3S)»  (39)  and  (40),  when  the  bang-bang  control  as  represented 

by  equation  (44)  la  utilized.  That  is,  the  control  function 

defined  by  equation  (44)  will  not  cause  the  representative 

point  to  remain  along  the  singular  control  surface  S  .   For  this 

s 


6b 


Fig. 28.  A  typical    optima!    trajectory  in      x,t    space 

showing     singular    and    nonsingular   subarcs  [9]. 
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reason,  It  is  easy  to  overlook  singular  solutions  when  ordinary 

computer  searching  methods  are  used  to  integrate  the  canonical 

equations  [12].  In  principle,  the  determination  of  the 

switching  surface  S  is  accomplished  by  expressing  the  initial 

value  of  the  (optimal)  adjoint  vector  z  as  a  function  of  the 

initial  conditions  vector  x(0).  In  practice,  this  procedure 

leads  to  a  rather  difficult  analytic  problem  which,  so  far,  has 

been  solved  in  only  a  few  special  cases.   Nevertheless,  If  the 

solution  for  the  optimal  control  involves  subarcs  of  singular 

control,  then  it  may  be  possible  (with  emphasis  on  the  word  "may") 

to  determine  z(x,  t)  on  S  without  knowing  the  general  expression 

for  z(x,  t),  by  making  use  of  the  fact  that  on  S  the  switching 

s 

function  and  all  Its  time  derivatives  are  identically  zero  [12]. 
However,  this  method  involves  the  heaviest  of  algebraic 
manipulations  even  for  two  dimensional  problems,  and  certainly 
encumbers  the  analytic  procedure  for  three  dimensional  problems 
such  as  the  present  case.  Also  the  fact  that  in  addition  to  the 
algebraic  manipulations,  a  certain  line  integral  has  to  be 
evaluated,  makes  the  above  mode  of  approach  to  determining  the 
singular  solution  quite  formidable. 

A  more  direct  numerical  approach,  as  outlined  by  Grethlein 
and  I.apidus  [15]  In  the  paper  entitled  "Time  optimal  control  of 
non  linear  systems  with  consti-aints",  was  used  to  obtain  the 
control  policy  for  the  present  problem.   It  should  be  borne 
firmly  In  mind  that  the  physical  realizabillty  of  singular 
control  situations  is  not  sufficient  to  establish  the  optircality 
of  singular  control  and  consequently  the  control  policy  as 
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obtained  by  the  method  of  Grethleln  and  Lapidus  may  not  be  a 
truely  optimal  policy  but  It  oertainly  is  one  of  the  methods 
available  for  determining  autoptical  policies. 

The  approach  taken  by  Grethleln  and  Lapidus  is  to  predict, 
over  one  sampling  period  by  means  of  the  na  their,  at  ical  process 
models,  the  responses  for  a  selected  number  of  control  levels. 
A  modified  objective  function  (usually  of  the  least  squares 
form)  Is  defined,  and  by  evaluating  this  objective  function 
for  the  predicted  values,  the  dptimum  control  action  is  selected 
at  each  sampling  period.  It  is  that  control  action  which 
minimizes  the  objective  function.  The  overall  control  system 
with  its  optimizing  scheme  is  called  the  "optimum  predictor- 
controller". 

In  general  the  dynamic  description  of  a  process  can  be 
represented  by  a  set  of  first  order  differential  equations 
written  in  vector  form  as 

^^[i(t)i  e(t)]  (M] 

where  x(t)  represents  the  state  variables  and  0(t)  represents 
the  control  variables.  Thus  the  state  of  the  process  is 
interpreted  as  a  vector  in  state  space  and  each  state  of  the 
process  which  Is  different  from  another  has  a  unique  set  of 
coordinates. 

In  principle  the  transient  behavior  is  obtained  by  solving 
the  differential  equations  simultaneously  with  the  proper 
initial  conditions  and  control  variables.   As  a  result  a  future 
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state  of  the  process  at  any  time  x  car.  be  evaluated  If  some 
initial  state  x(t.)  is  given  in  addition  to  the  control  vector 
0(t)  for  all  time  t  <  t  <  T.  The  existence  of  a  solution  does 
not  imply  that  it  is  necessarily  an  analytic  solution.  In 
most  nonlinear  cases  only  numerical  solutions  are  obtainable. 

Since  a  computer  can  only  operate  on  digits  vihich  represent 
the  state  of  the  process  at  a  discrete  point  In  time,  any  control 
system  utilizing  a  digital  computer  will  necessarily  be  a 
sampled  data  system.  From  the  point  of  view  of  the  controller, 
a  process  Is  represented  as  a  sequence  of  numbers  spaced  in  time. 
Similarly  the  control  action  is  a  sequence  of  numbers.  It  is 
convenient  to  take  the  time  Interval  between  sampling  points 
equal  to  X  so  that  real  time  can  be  represented  at  the  sampling 
points  by 

t  -n   X ,  2\ ,  3X ,  . . . ,  kX  ,  ... 

Althoiigh  the  control  and  state  vectors  of  the  process  vary 
continuously  with  time,  their  specific  value  at  the  sampling 
point  t  =  kX  is  represented  by  6(kX)  and  x(kX)  respectively  or 
more  simply  by  6(k)  and  x(k).  If  the  sampling  period  is  taken 
small  enough  so  that  the  control  0(t)  for  kX  <  t  <   (k  -i  1)X  can 
be  considered  essentially  constant  and  equal  to  e(k)  for  the 
entire  period,  the  dynamic  equation  for  the  process  will  become 

£§=<£[>;  6(k)]  (jty) 

for 

kX  <  t  <  (k  +  1)X. 
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Thus  the  dynamic  behavior  of  the  process  can  be  expressed  in 
the  form 

x(k  +  1)  a£[x(k)|  o(k)j  X]  (40) 

With  a  knowledge  of  the  state  of  the  process  and  the  input  vector 
at  any  time  t  a   k,  equation  (43)  gives  the  state  of  the  process 
at  the  next  sampling  point. 

Before  it  is  possible  to  single  out  a  particular  controller 
as  being  optimum,  it  is  necessary  to  define  quantitatively  some 
measure  of  merit  by  which  the  performance  of  the  system  is 
evaluated.  The  simplest  type  of  measure  of  performance  is  some 
function  of  the  difference  between  the  actual  output  of  a  process 
at  a  given  time  and  that  desired  for  the  dynamic  system.   A 
convenient  measure  may  be  defined  as 

S(t)  *  a1[xi(t)  -  a*]  +  a O  (t)  -  x*]  +  ...      (49) 

where  a  ,  a  ,  ...  are  appropriate  weighing  factors  and  x..  ,  x*,  .,, 

are  the  coordinates  of  the  desired  final  state.  Assuming  that 

the  control  vector  is  changed  only  at  the  sampling  point, 

equation  (49)  is  summed  over  all  time 

N 
J(N)  a   £   S(kX)X        N  — >  oo  (50) 

k=l 

where  N  is  taken  sufficiently  large  to  cover  the  entire  transient 
period.  Thus  the  optimum  dynamic  response  is  obtained  when  a 
discrete  set  of  control  inputs  e(k)  are  found  such  that 
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N 
J  (H)  =  Kin   I  8(kX)       N  --—>«>  (51) 

6(k)  k=l 

At  this  stage  mention  must  be  made  that  if  the  system  is  table, 

the  dynamic  response  will  ultimately  come  to  an  equilibrium  or 

steady  state  at  some  k  -■.   NeqUii" 

The  optimum  control  action  at  any  given  sampling  period  is 

generated  by  the  following  computation  scheme.  Since  the  present 

state  of  the  process  x(k)  is  known  as  a  result  of  a  feedback 

measurement ■  the  state  of  the  process  at  the  next  sampling 

period  x(k  •;•  1)  can  be  calculated  from  equation  (^8)  for  any 

number  of  control  actions.  Specifically,  the  state  of  the  process 

Is  calculated  for  the  maximum  control  G   ,  the  minimum  control 

max 

action  0   ,  and  three  intermediate  control  actions  6,,  8  and 
niin  J-   <!■ 

6  ;  thus 

3  ' 

,   xa(k  +  i)  =^[x(k),  emaxi  x] 

x2(k  +  1)  -_-<£[x(k)j  8   i  X] 
T  L      min 

x3(k  +  1)  =  <£[>(k);  8  |  X]  (32) 

*(X   +  1)  =  £|>(k)i  62s  X] 

x5(k  +  1)  =^[x(k);  Go!  X] 

where  the  superscript  numbers  are  used  to  distinguish  the 
predicted  states.   The  responses  to  the  various  control  actions 
are  evaluated  with  the  aid  of  the  objective  function  for  the 
(k  +  1)  period,  namely 

P1(k  +  1)  =  P[x*i  xX(k  +  l)j  (53) 
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P'(k  +  1)  =.  P[x*;  S*(k  +  1)] 

P3(k  +  1)  =  P[x*j  x3(k  H  1)] 

(53  cor.t.) 
P  (k  +  1)  =  P[x*j  x^k  +  1)] 

P5(k  +  1)  =  P[x*;  x5(k  +  X)] 

The  possible  configurations  for  the  predicted  objective  function 

and  the  control  action  are  shown  in  Figure  29.  The  computer 

routine  selects  the  optimum  control  action  at  t  =  kX  on  the 

basis  of  a  predicted  minimum  objective  functions  at  t  a   (k  +  1)X. 

For  configuration  (a)  or  (b),  9    is  selected  as  the  optimum 

max 

control  action.   For  configuration  (c)  or  (d),  6    is  selected 

min 

as  the  optimum  control  action.  The  need  to  predict  the  performance 
for  some  intermediate  control  action  becomes  clear  when  the 
optimum  is  not  one  of  the  extreme  values.  When  configuration 
(e),  (f)  or  (g)  occurs,  the  computer  fits  a  second  order  curve 
through  the  minimum  point  and  its  two  neighboring  points  already 
calculated,  and  finds  the  minimum  in  the  curve.  The  corresponding 
value  of  9  becomes  the  optimum  control  action  for  that  period. 

The  computer  logic  flow  chart  for  the  above  method  is 
shown  in  Figure  30.  The  explanation  of  the  various  steps  in  the 
corresponding  numbered  boxes  in  Figure  30  is  as  folloifs. 

1.  The  five  values  of  the  decision  variable  are  read  in, 

where  AL(1)  =  L  .   AL(5)  -.-  L    and  AL(2),  AX.(3),  and  AL(4) 

min         max 

are  three  intermediate  values  of  the  decision  variable  L, 

?..   The  lower  boun;l  of  the  first  sampling  period  (PRNT(U), 
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Fig. 29.    Possible      configurations    for    predicted 
performance     function. 


72 


/2T\  M 

<DEV(J)<DMi^>J     JMIN=  J 

\  y  « — — ' 


Fig.30.     Computer    logic    flow    chert   for    method   of 
Grefhlein    ond  Lopidus. 
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D|=(R|-R2)  (RI-R3) 
D2=(R2-R3)  (R2-RI) 
D3*(R3-R)      (R3-R2) 
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A  =  EI+E2+E3 

B»EhGI  +  E2    G2  +  E3  G3 

C  =  E1+HI  +  E2    H2+E3    H3 
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the  upper  bound  of  the  first  sampling  period  (PHMT(2)),  and 
the  integrating  increment  (PBMT(3))  are  initialized. 

3.  The  number  of  sampling  periods  (KN).  the  current 
sampling  period  (N) ,  and  the  Initial  minimum  deviation  (SUMT) 
are  initialized.  The  switch  II  is  set  equal  to  unity.   (II. 
can  be  either  1  or  0). 

4<   Tlie  subscript  J  is  Initialized. 

5.  The  decision  variable  (AP)  is  initialized. 

6-8.   The  current  value  of  the  sampling  period  is  checked 
and  accordingly  the  initial  values  of  each  sampling  period  are 
specified.  For  N  n  1,  the  Initial  values  are  specified  in  Box 
8,   For  N  >  1,  the  initial  values  are  specified  in  box  7  wherein 
the  terminated  state  of  the  Nth  sampling  period  becomes  the 
initial  state  of  the  (N-j-l)th  sampling  period. 

9.  Subroutine  RKGS  is  called  and  the  integration  commences. 

10,  11.  The  extent  of  the  integration  is  checked  to  see 
whether  the  upper  bound  of  the  current  sampling  period  has  been 
reached  if  so,  command  transfers  to  box  12,  and  if  not,  command 
passes  to  box  11  where  the  time  is  incremented  and  the 
integration  proceeds. 

12.  The  value  of  the  switch  II  is  checked.   If  II  =  1, 
command  transfers  to  box  13;  otherwise  to  box  41, 

13.  The  quantity  SUM  is  defined  as  the  minimum  deviation 
up  to  the  (K-l).th  sampling  period  (3UHT)  plus  the  current 
deviation  (S) . 

14.  The  variable  DEV(J)  is  set  equal  to  SUM  and  this 
represents  the  total  deviation  corresponding  to  the  curz-ent  value 
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of  J. 

15.  -  DEV(J)  is  written  out. 

16,  17.  The  subscript  J  is  checked  to  make  sure  all  the 
5  decision  variables  defined  in  box  1  have  been  used.  If  J  is 
less  than  5,  it  is  incremented  by  unity  and  command  is  passed 
to  box  5.   If  J  is  equal  to  5  command  passes  to  box  18. 

18-24.  The  minimum  of  the  five  DEV(J)  is  determined,  and 
the  value  of  J  corresponding  to  the  minimum  DKV(J)  is  stored  as 
JMIN. 

25.  26.   The  new  variables  AL1  to  AL5  and  DEVI  to  DEV5 
are  defined.. 

27-30.   If  JMIN  =  1,  then  ALEM  is  defined  as  the  decision 
variable  corresponding  to  JMIN  and  DEVJ  is  defined  as  the 
deviation  corresponding  to  JMIN.  DEVJ  is  written  out  and 
command  is  transferred  to  box  38. 

31-34.  If  JMIN  =  5,  then  ALEM  Is  defined  as  the  decision 
variable  corresponding  to  JMIN  and  DEVJ  is  defined  as  the 
deviation  corresponding  to  JMIN,  DEVJ  is  written  out  and  command 
is  transferred  to  box  38. 

35-3?.  If  JMIN  =  2,  3,  or  4,  command  is  transferred  to 
box  44,  45  or  46  respectively  whereby  a  second  order  curve  will 
be  fitted  as  shown  in  Figure  29  (e),  (f)  and  (gj  and  the  minimum 
of  the  curve  determined. 

38.  The  decision  variable  is  set  equal  to  the  minimizing 
decision  variable  as  determined  from  boxes  28,  29,  or  50. 

39.  SUKT  is  set  equal  to  DEVJ  as  the  minimum  total 
deviation  up  to  the  present  sampling  period. 


7? 

40.  The  switch  II  Is  set  eqml  to  zero  and  command  Is 
transferred  to  box  6.   The  reason  for  setting  II  =  0  is  that 
once  the  minimizing  decision  variable  has  been  determined  then 
the  computer  can  skip  all  the  command  from  boxes  13  to  40 
inclusive,  thereby  command  being  transferred  from  box  12  to  box 
41. 

41.  The  values  of  the  state  variables  x.  ,  x  and  x  are 

12      3 

written  out  along  with  the  minimizing  decision  variable  end  the 
number  of  the  current  sampling  period. 

42.  Herein  are  initialized  the  various  parameters  for  the 
next  sampling  period  and  II  is  again  set  equal  to  1. 

43.  If  the  number  of  the  new  sampling  period  has  not 
exceeded  the  value  of  NN  then  command  is  transferred  to  box  4 
and  the  above  procedure  repeated.  If  the  number  of  the  new 
sampling  period  exceeds  NN,  the  computation  in  terminated. 

44-49.  The  coefficients  of  a  second  order  polynomial 

through  three  points  are  calculated.  The  polynomial  is  of  the 

2 

form  y  =  ax  -  bx  +  c  where  y  corresponds  to  the  deviation  and 

x  corresponds  to  the  decision  variable. 

50.  The  minimizing  decision  variable  is  determined  by 

setting  d£  _  0,  I.e.  2ax  -  b  =  0,   x  =  b/2a. 
dx 

51.  From  the  minimizing  decision  variable,  the  corresponding 
minimum  deviation  Is  determined. 

52.  The  minimum  deviation  DEVJ  is  written  out  and  command 
is  transferred  to  box  38. 

The  above  method  was  applied  to  Cases  1,  2,  3  and  4  and  the 
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optimal  (or  suboptimal)  policy  was  obtained  In  each  case.  In 

order  that  the  method  be  applicable  to  our  problem  a  few  changes 

were  incorporated!  namely  the  control  action  was  selected  so  that 

the  dynamic  response  would  minimize  an  objective  function  of  the 

.type  given  in  equation  (49)  rathers  than  the  original  form  in 

equation  (33).  Since  in  the  present  example  concentrations  are 

being  controlled  a  suitable  overall  objective  function  would  be 

that  defined  by 

N 
J(N)  =  £  S(k)X 
k=l 

with 

2  2  2 

S(k)  »  aCx1(k)  -  x*]  +  b[x2(k)  -  x*]  +  eOgdO  -  x*] 

(5*) 

where  x.  x  ,  and  x  are  the  coordinates  of  the  desired  final 
12      3 

state  S'  5  (x  =  0.71399,  x  =  0.60966,  x  =  0.2824?) ,  and  a,  b, 

and  c  are  appropriate  weighting  factors.  According  to 

Grethlein  and  Lapidus  [1 5]  there  is  no  prior  way  of  determining 

the  exact  weights,  but  the  dynamic  properties  of  the  system  can 

serve  as  a  guide  line.  In  our  particular  problem  the  system  is 

to  be  transferred  from  an  Initial  state  S*  "=  (x  =  0.7146?, 

x  =  0.60909,  x  =  0.24899)  to  the  terminal  state  S'  e  (x  = 

0.71399.  x  a   O.60986,  x  -  0.2824?)  in  a  time  optimal  manner. 

Comparing  S'  and  S'  we  see  that  x^  and  x„  don't  have  as  much 
12  12 

distance  to  cover  in  the  state  space  as  x  does  and  consequently 
we  weigh  the  deviation  of  x  more  heavily  that  the  deviation  of 
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x  and  x  .  The  rationale  for  doing  this  Is  that  in  view  of  the 
larger  margin  of  operation  in  x  and  hence  a  larger  margin  for 
its  deviation,  we  penalize  it  heavily  thus  making  sure  that  it 
attains  the  desired  final  state  without  undue  deviation.  Thus 
the  control  action  is  biased  mainly  towards  the  deviation  of 
x  .  In  making  x  attain  its  desired  final  state  in  the  most 
direct  manner  we  also  ensure  that  x  and  x  attain  their  desired 
final  states  since  x  ,  x  ,  and  x  are  interrelated  by  material 

J.  c  } 

balances  and  also  the  system  dynamics  are  based  on  these 
material  balances. 

The  predictor  controller  that  is  used  along  with  an  objective 
function  of  the  type  in  equation  (5*1),  makes  certain  that  the 
system  attains  the  final  state  S'  for  a  wide  range  of  values  of 
the  weighting  factors  a,  b  and  c.  However  since  we're  Interested 
in  time  optimality  we  must  select  a,  b  and  c  such  that  S*  is 
a^'.alned  in  the  shortest  possible  time. 

In  experimenting  with  the  above  controller  it  was  found 
that  the  response  time  T  in  going  from  S*  to  S'  was  significantly 
decreased  as  the  ratio  of  cia  and  cib  was  increased.  Thus  the 
weighting  factors  were  picked  such  that  a  and  b  are  equal  and 
the  ratio  cia  or  cib  is  an  integral  power  of  10. 

The  scheme  used  for  selecting  the  optimal  values  of  a,  b 
and  c  is  as  follows i 

1.  Given  the  initial  state  S'  the  optinal  predictor 
controller  was  utilized  so  as  to  minimize  the  objective  function 
in  equation  (5*0  with  ratio  of  weighting  factors  cia  -■  cib  -■   10, 

2.  The  time  at  which  the  desired  state  S'  is  attained 

2 
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is  noted  as  T,  « 

3.   The  procedure  In  step  1  is  repeated  with  eta  =  ctb 
=  10  for  1  =  2,  3,  ...  and  the  time  at  which  state  S'  3k 
attained  is  noted  as  f  ,  ,  i  =•  2,  3i  ••• 

The  above  listing  of  X      gives  the  times  required  to  transfer 

the  system  form  state  S*  to  the  state  S'  for  increasing  values 

J.  *- 

of  the  ratios  eta  or  cib.   It  was  found  that  as  i  increases  X 
decreases  and  that  there  is  a  certain  1  beyond  which  any  further 
increase  in  1  does  not  produce  any  decrease  in  f  .   This  is  the 
1  that  determines  the  optimal  ratio  eta  =  ctb  B-  10  .   These 
values  of  a,  b  and  c  and  the  corresponding  control  policy  which 
transfers  the  system  from  S'  to  S'  was  then  accepted  as  the 
optimum  predictor  controller  achieving  time  optimality.   The 
value  of  1  and  hence  the  values  of  a,  b  and  c  are  found  to  vary 
as  we  go  from  Case  1  to  Case  k. 

Hence  in  recapitulating  we  must  realize  that  the  original 
objective  function  in  equation  (33)  which  demands  time  optimality 
has  been  discarded  and  in  its  place  a  modified  objective  function 
of  the  type  in  equation  (^9)  is  adopted.   The  weighting  factors 
in  equation  (49)  are  then  experimented  with  so  that  the  system 
is  transferred  from  the  state  S'  to  the  state  S'  in  the  most 
direct  way  8nd  In  a  time  optimal  manner. 

The  optimal  response  for  the  four  cases  is  shown  In  Figures 
31.  32,  33,  and  3^»   Basically  two  types  of  control  policies 
are  encountered.   The  first  type  (in  cases  1,  2  and  3)  wherein 

initially  L    is  used  until  the  first  switch  after  which 

max 

singular  control  is  used  finally  ending  up  asymptotically  at  L'. 
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The  second  type  (in  Case  k)   uses  L  _      initially  until  the  first 

switch  after  which  L  .   is  used  until  a  second  switch  after 
min 

which  singular  control  is  used.   In  all  four  cases  the  control 
variable  L  tends  to  a  final  value  of  L'  and  as  a  result  the 
system  approaches  the  desired  final  state  saymptotieally  without 
overshoot  (in  contrast  to  the  case  of  the  reference  system  where 
purely  bang  bang  control  was  used)  and  maintains  the  desired 
final  state.  The  results  for  the  four  cases  are  tabulated,  in 
Table  2.   In  Figure  35  a  comparison  is  made  between  the  bang- 
bang  policy  as  obtained  from  the  reference  system  analysis  and 
the  singular  policy  as  obtained  from  the  mixing  pool  model 
system  analysis,  for  the  time  optimal  problem. 

The  method  of  Grethlein  arid  Lapidus  [15]  was  also  applied 
to  the  time  optimal  problem  for  the  Reference  system  treated  in 

Chapter  2.   The  original  objective  function  In  equation  (12), 

T 
I.e.  /  dt,  was  replaced  by  a  weighted  least  squares  objective 

o 
function  and  the  procedure  described  above  was  carried  out  to 

obtain  a  tine  optimal  policy.  These  results  are  compared,  in 
Table  3,  with  the  rigorously  determined  bang-bang  policy  of 
Chapter  2,  end  the  corresponding  responses  are  compared  in 
Figures  36,  37,  38  and  39.   It  is  quite  evident  that  the  bang- 
bang  policy  does  indeed  achieve  time  optimality.  The  responses 
from  the  two  policies  differ  only  in  that  the  singular  response 
approaches  the  final  desired  state  asymptotically,  while  the 
bang-bang  response  tends  to  overshoot  at  the  final  state. 
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TABLE  3 

COMPARISON  OF  RESULTS  FROM  BANG-BANC,  (OBJECTIVE 

T 
FUNCTION  IS  /  at)  AND  SINGULAR.  (OBJECTIVE  FUNCTION 

o 
IS  OF  THE  WEIGHTED  LEAST  SQUARES  TYPE)  POLICIES  OF 
THE  TIME  OPTIMAL  PROBLEM  FOB  THE  REFERENCE  SYSTEM. 


Bans  Bang 

S  •  ■,:<;<•''   T.i' 

Final  time 
T,  mln 

Final  time 
T ,  mln 

Case  1 

2.8 

4.85 

Case  2 

4.05 

5.45 

Case  3 

8.16 

8.46 

Case  4 

1.68 

4.3 
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Fig.  35.   Comparison    of   bang  bang  (reference  system)  and 
singular  (two  tanks  in  series  system)  control  policies 
for  the    time   optimal   problem. 
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CHAPTER  ■': 


CONTROL  WHICH  IS  OPTIMUM  IN  THE  SENSE  OF  MINIMUM  TOTAL 
DEVIATION 


In  this  chapter  attention  Hill  be  fooussed  on  the  minimi: 

tion  of  an  objective  function  of  the  form 

T  ,        2 
9  «=  /  [x*  -  x(t)]  at 

o 

wjilch  is  a  measure  of  the  deviation  from  some  desired  quantity 
x  over  the  unspecified  time  interval  T,  A  control  policy  L(t) 
ie  sought,  which  transfers  the  system  from  some  known  initial 
state  to  a  desired  final  state  such  that  S  is  minimised. 

The  two  systems  which  will  be  subjected  to  such  a 
minimization  are  the  reference  system  and  the  system  with  the 
mixing  condition  on  the  top  tray  described  by  the  two-tanka  in- 
s cries  model. 

1.   REFERENCE  SYSTEM 

In  this  section  we  shall  consider  the  reference  system 
only.  The  performance  equations  of  the  system  are  those  of 
Chapter  2  and  are  repeated  here  for  convenience 

dt  "      KT      XX   +  HTX2  +  HfJ,  «« 

B  B  B 

The  boundary,  conditions  are 


y... 


x  (0)  =  O.63OO        x  (TJ  =-  0.6300 

x  (0)  =  0.2767        x  (T)  =  0.3067 
2  2 


where  T  is  unfixed. 

The  problem  is  to  determine  the  control  L(t)  such  that  the 
system  is  transferred  from  the  initial  state  S  =  (O.63OO, 
0.2767}  to  the  final  state  S  a  (O.630O,  O.3067)  simultaneously 
minimising  the  objective  function 

m 

S  ,-  /  [z*  -  x  (t)]2  dt 

where 

x*  =,   0.837  [cf.    equation   (7)] 

x  (t)  =  y    =-•  0.4b  (t)  +  0.56  [cf.    equation   (7)] 

t 
and  the  final  time  T  is  unspecified.   Thus  we  essentially  want 

to  minimize 

T 

0' 

Introducing  an  additional  state  variable,   wa  have 


S  ,=     /      [0.837  -  O.Wx,(t)    -   0.56]?~  dt  (57) 

n  1 


T  2 

x   (T)   •.=   S'   -_■:      /    [0.837   --  0.'44x    ft)    -   O.56]   dt 
3  0  1 


dx,  2 

-%£  =   [0.837   -   O.^   -   0.56]    ;      x3(0)   =   0  (58) 


The  Haralltonlan  for  the  system  becomes 
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„   /Lra  -  Vm  -  Li-.  „   .  Vm. .  „    Lc     L 
T  T        T      B 

/ Vm  +X tJL." _ir'  ';  — _-  -r 

B  B 

+  [0.63?  -  m^  •-  o]  a3  (59) 

The  adjoint  equations  can  be  derived  from  the  Hamiltonian 

<g  m  .   (^  -^  '  ^  -  ^  ♦  2^(0.837  -  Ml  -  o)Z3 


(60) 


dz 

IF  = 


B  -  +  (Vm  +FHtL-  Y}  Zg  (6l) 

T  B 


dz, 

-jg   -  0  (62) 


The  boundary  conditions  on  the  adjoint  variables  are 


z  (0)  unspecified  Z (T)  unspecified 
z  (0)  unspecified  z„(T)  unspecified 
8,(0)  unspecified       z  (T)  =  0 


Equation  (62)  along  with  the  boundary  condition  on  z  Implies 
that 


z  «=  1,      0  <  t  <  I 

Making  use  of  the  property  that  the  minimum  of  the  Hamiltonian 
Is  zero  when  T  Is  unspecified,  we  have 


mzlXl       x1z1       Czx       Xl?2       x2z2  Vm  xgZi        Vm  xlZ% 

HT  HT  h,r,  HB  Hfi  HT  hT 


Vm  x2Z2        (F  -   V)X222        (Fxp  -   Vo)Zg 

.......     ...   _~      :    +  _.      _ 

BE  B 


+  [0.837  -  mx  -  c]2  „  0  (63) 


wherein  the  Hamiltonian  has  been  rearranged  and  shown  linear 
in  the  control  variable  L,   Once  more  a  bang-bang  policy  is 
indicated  vihereby 

L  a  L    if  —-L-1-  -  -3L1  +  — 1  4  -1—2-  -  -?~-2)  <  0     (6ftJ 

max      h       H      H      II      H 
flT       T      T      B     HB 

or 

mln      HT      ET     HT    Hr     Hr 

Comparing  equations  (63)  and  {6k)   with  equations  (21)  and  (21a) 

respectively  it  appears  that  the  minimum  total  deviation  and 

time  optimal  problems  have  similar  solutions.   As  a  matter  of 

fact,  the  performance  equations  being  the  same  for  the  two 

oases,  the  phase  planes  are  identicaland  hence  the  bang-bang 

policy  obtained  for  the  time  optimal  problem  seems  to  be  feasible 

for  the  minimum  total  deviation  problem.  However,  this  is  not 

so.  The  trial  and  error  method  outlined  in  Figure  3  was  applied 

to  the  present  problem  wherein  the  initial  value  of  z  was 

2 

assumed  and  the  initial  value  of  a  was  obtained  froa 

z   (0)   -   -       0.HH3  L(0)   -  0.1213  ,    ■  C65) 

VU;   ~  -0.2071   L(o)    -   0.2071      V0'  V    5' 


9? 

Equation  (65)  is  the  analog  of   equation.  (22)  in  the  time  optimal 
case.   A  bang-bang  polio/  vias  obtained  for  the  four  cases 
investigated  in  Chapter  2  i.e. 


Case  1   L  ,   =  0.833     L    =  0.997^ 
min  max 


Case  2   L  ,  =  0. 86588   L    =  0.96'l-52 
min  .  max 


Case  3   L  ,     *-.   0. 89876   L    =  0«93l6^ 
min  max 


Case  I*       L  .   =  O.833 
min 


I.    =  1.333 

max 


however,  the  Hamlltonian  retained  a  positive  value  along  the 
optimal  response,  instead  of  attaining  a  minimum  at  zero  for  all 
four  cases,  thus  indicating  that  the  objective  function  could 
not  be  minimized  by  mere  bang-bang  control.   Once  more  a  singular 
policy  was  obtained  by  using  the  method  of  Grethlein  and  Lapidus 

CV]. 

In  order  to  apply  the  singular  control  method  a  modified 
objective  function  was  defined  of  the  form 


N 
J(N)  -_-  Z   S0O 

toal 


H  ■— Js>  °° 


(66) 


with 


S(k)  a   [0.83?  --  O.kk  Y      -  0.56]' 


(6?) 


Equation  (67)  should  be  Interpreted  as 

,2 


s(k)  -,  [x:  -  x  (t)] 


where 
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x  =  0.837 
D 


and 


x  (t)  =  O.k'i   x  -  O.56      [cf.  equation  (?)] 


Note  that  in  equation  (6?)  the  weighting  factor  is  taken  as  unity 
The  sequence  of  control  variables  vias  selected  so  as  to 
minimize  equation  (6?)  at  each  sampling  point  in  tine.  Also, 
the  variable  x  (t)  as  defined  by  equation  (58)  was  evaluated  all 
along  the  optimal  trajectory  (optimal  with  respect  to  equation 
(66)).  Becall  that  x  (T)  (T  is  the  time  at  which  the  final 
desired  state  is  reached)  is  indeed  the  objective  function  that 
we  originally  wished  to  minimize  according  to  equation  (58). 
It  vias  found  that  x  (T)  as  obtained  by  the  singular  policy 
is  indeed  lower  then  Xo(T)  as  obtained  by  the  bang-bang  policy. 
Comparison  of  the  two  is  made  in  Table  k,   for  all  four  cases. 
The  response  to  the  singular  policy  is  shown  in  Figures  ho,   k\, 
kZ   and  43.   The  response  to  the  bang-bang  policy  is  identical 
to  that  In  Figures  15.  16,  17  and  18.  In  Figure  kk,   comparison 
is  made  of  the  control  policies  resulting  from  bang-bang  and 
singular  control.   Rather  significant  differences  can  be  seen  in 
the  form  of  control  as  well  as  the  duration  of  control. 
Though  acheivlng  a  smaller  deviation,  singular  control  must  be 
applied  for  a  longer  time  while  bang-bang  control  which  does 
not  acheive  as  small  a  deviation  as  singular  control  is  applied 
for  a  very  short  time.   It  seems  as  though  in  a  practical 
situation  the  most  likely  objective  function  to  minimize  would 


99 


TABLE  4 


COMPARISON    OF   RESULTS    FROM   BANG-BANG   AND  SINGULAR  POLICIES 
FOR  MINIMUM   DEVIATION   PROBLEM    (REFERENCE   SYSTEM) 


Bang-Bang  Singular 

Total  minimum     Final   time  Total  minimum     Final   time 

deviation             T,   ruin.  deviation             T,   rain 

x3(T)  x   (T) 


-4  -6 

Case  1     0.2154  x  10  2.80  O.38OO   x  10  ?.?0 

Case  2     0.1034  x  10  4.05  O.6730   x  10"  7,80 

-4  -4 

Case  3     0.2583  x  10  8. 16  0.248?1  x  10  9-90 

Case  4     0.1031  x  lo"3  1.68  O.368O  x  10"6  7. 60 


a  oo 


Time,  minutes 

Fig. 40.     Response    of    reference     system     to    the 

singular    policy     for    the     minimum     deviation 


problem  (Case  I). 
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0-640 

x'  0  ■  6  3  0 

0-620 

0-3  !  0 
0-300 
0-290 
0  280 


x 


0-96452 


-J   0-9152 


0  86  588- 


_L_ 
2 


4  e  8 

Time,     minutes 

Fig. 41.    Response     of    reference    system     to    the    singular 
policy    for    the    minimum    deviation    problem  (Case 2) 
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0-93164 
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-J     0-91520 

— 

^ 

0-89876 

„  J 
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■           i 

2 

4 

6                  8 

Time,    minutes 

Fig. 42.    Response    of    reference    system    to   the    singular 

policy    for    the    minimum     deviation    problem(Case3). 
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4  6 

Time,      minutes 


Fig. 43.    Response    of    reference     system    to    the 

singular    policy    for   the    minimum     deviation 
problem  (Cose  4). 
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0-9974 


0-915?. 
0-8330 

4  6 

time,minuies. 

Fig  .44  Comparison   of  bang  bang    end  singular  control   policies 
os   applied    to    the    minimum    deviation   probiem    and 

tine    reference   system. 
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be  of  the  type 

(t)]~  +  1\   dt 


slH  -  xM  *  J 


wherein  a  compromise  Is  made  between  minimum  time  and  minimum 
deviation. 

2.   SYSTEM  WITH  TOP  TRAY  DESCRIBED  AS  TWO  TANKS  IN  SERIES 

The  performance  equations  of  Chapter  3  are  repeated  here 
for  convenience. 

uxl    ,Lm  -  Vm  -  2Lv,   .  Ltn   .  Vm     2Lc         ,-„, 

Tt    B  < 5T~ )X1  +  H~X2  *   H"x3  +  H~         (68) 

T  T      T      T 

dx 

2   2L  _    r2L  +  Vm>    .  Vm  fiCa, 

__  .  --  x     -   ( *—~ )x     +  —  x  (69) 

T  T  T 

dt    -   H  x2        l  H  ;    3  +  H  uu' 

B  B  B 

c 

and  the  boundary  conditions  are 

x  (0)  B  0.71W7  x  (T)  -.,  0.71399 

x   (0)  =  0.60909  x  (T)  =  0.60986 

X   (0)   *   0.2^99  x  (T)  =  0.232^7 


where  T  is  unfixed. 

The  problem,  once  again,  is  to  determine  the  control  L(t) 

which  will  transfer  the  system  from  the  initial  state  S'  = 

1 
(0.71^67,  0.60909,  0.2W399)  to  the  final  state  S'   5   (0.71399, 

0.60986,  0.282'J-7)  and  in  so  doing  will  minimize  the  objective 
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function 


S  =  /  [x*  -  r  (t)]  dt 


where 


x*  =  0.85123       [of.  equation  (2?)] 

x  (t)  »  0.22x  +  0.22x  +  0.56       [cf.  equation  (27)] 

D  A.  C. 

and  the  final  time  T  is  unfixed.  Thus,  essentially  we  have  to 
minimize 

S'  =  /  [0.85123  -  0.22X,  -  0.22x„  -  0.56]2  dt       (71) 
0  12 


Next  a  new  variable  x.  is  introduced  such  that 

T 

x. (T)  „  S'  m     !     [0.85123   -  0.22X.    -   0.22xo  -   0.56]Zdt 
4  0  12 

or 

dxa  2 

-d-~  =  [0.85123  -  0.22x2  -  0.22x2  -  O.56]  j  xk(0)   =  0 


(?2) 


The  Hamiltonian  for  the  system  becomes 


H  .  (£"  -  Vm  -  2L)Xi2^  imXoZi+  gn^^  2Lc 
'  T 


"IT    H      2T  HnTl     H        1 


■*■  J-  J  H 

/Vm  +F+   L  -   V\  ,FxF  "   Vc, 

-   (._ ±___ )xf3  +    (      '^      )    ,3 

2 
+   [O.85123  -   0.22xn    -   0.22x„   -   0.56"J    z.  (73) 
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and  the  adjoint  differential   system  derived  from  the  Bafiiiltonian 
is, 


dzl  fLm  -  Vm  -  2LW 


+  O.W[0. 85123  -  0.?2x1~  0.22x2~-  0.56]z^  (7';-) 

dt    -        H^l  +    l      HT        ;Z2        HB  3 

4    0.44[ 0.85123   -   0.22x     -   0.22x     ■•   0.56]z  (?5) 


dt    -        H  *1       H  z2  +    K  H  ;23  uo; 

T  T  B 

-dF  -  °  CWJ 

and  the  boundary  conditions  on  the  adjoint  variables  are 

z  (0)  unspecified  z   (T)  unspecified 

z  (0)  unspecified  z„(T)  unspecified 

z  (0)  unspecified  z  (T)  unspecified 

z,  (0)  unspecified  z,  (T)  „  1 

From  the  final  condition  on  z  and  equation  (?4),  it  follows  that 

z  =  1,     0  <  t  <  T 
The  Hamiltonian  can  be  rewritten  as  follows, 

H  =  (mx1z1  --  2x1^  +  mrgZi  +  2cz.x  +  2x  zg  -  2x  z 


(73) 
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jj.   _^..^i  _  _-^.-^}l  +  V-Tijc.z.-  Vmx-z.-  Vmx0z„+  Vmx„  z,. 


+  L0. 85123  -  0.22x  ••  0.22x  -  O.56] 


wherein  the  H_  has  been  omitted  since  it  has  a  value  of  unity,. 

Also  the  Hamiltonian  is  shown  linear  in  the  control  variable 

L,  thus  implying  a  bang-bang  policy.  The  switching  function  of 

the  bang-bang  policy  Is  Identical  with  equation  (kk)   which  is 

not  unexpected.  From  the  analysis  of  phase  planes  in  Chapter  3 

it  was  shown  that  the  final  state  could  not  be  reached  by  purely 

bang-bang  control  and  hence  in  the  present  problem  we  straightaway 

resort  to  singular  control. 

In  order  to  apply  the  singular  control  method  a  modified 

objective  function  was  defined  of  the  form 

N 
J(N)  =  Z   S(k)  N  — &■  ~  (80) 

k=l 

with 

S(k)  =  [0.85123  -  0.22Xl  -  0.22x2  -  0.56]2  (81) 

Equation  (81)  should  be  interpreted  as 

8(k)  ,  [>*  -  xD(t)]2 


where 


and 


xd  =  0.85123 
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X  (t)  =  0.22x  +  0.22x  +  0.56     [of.  equation  (72)] 

Note  that  in  equation  (81)  the  weighting  factor  has  been  taken 
as  unity.   The  sequence  of  control  variables  was  determined 
so  as  to  minimize  equation  (81)  at  each  sampling  point  In  time. 
Also,  the  variable  x, (t)  as  defined  by  equation  (72)  was 
evaluated  all  along  the  optimal  trajectory  (optimal  with  respect 
to  equation  (80)).   The  following  four  cases  were  investigated 

Case  1.   L-    =0.833     L=  1.011 
min  max 

Case  2.   L  .  a  0.8686    L    »  0.9752* 
min  max 

Case  "3.   Lmln,0.90to    L^  .  0.9398 
Case  4.   L^-  0.833     Lffiax  -  1.333 

and  the  corresponding  values  of  x. (T)  are  presented  in  Table  $, 
The  response  of  the  system  for  the  four  cases  is  plotted  la 
Figures  4-5,  46,  47  and  48. 
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TABLE 


RESULTS  OF  MINIMUM  DEVIATION  PROBLEM  FOR  SYSTEM 
WITH  TOP  TRAY  DESCRIBED  AS  TWO  TANKS  IN  SERIES 


Final  time 

Minimum  deviation 

T,   mln. 

x4(T) 

Case  1 

9-50 

0.1   x  10-8 

Case  2 

9.60 

0.227  x  10 

Case  3 

10.95 

0.2.5   x  10"'4, 

Case  k 

9.40 

0.1   x  10"9 

Ill 


! 

2 


4 


8 


time,  minutes 

Fig. 45.      Response       of      the      two     tonlts      in      series      system 
to     singular      control     for     the      minimum    deviation 
problem    (Case.  !.) 
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0-2S0 
0-270 
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6 


2  4  6  8  10 

Time,     mirv.ites 

Fig. 46.    Response    of    the    two     tanks     in    series 

system     to     singular    control    for    the     mini- 
mum    deviation     problem  (Case  2  ). 
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Fig. 47.     Response    of    the     two     tanks     in    series     sys- 
tem    to    sincjulGr    contro!    for    the     minimum 
deviation     problem  ( Case ■  3.) . 
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Fig. 48.      Response      of    the    two     tanks    in 
series     system     to     singular     cont- 
rol    for    the     minimum     total      de- 
viation    problem. (case  4.) 
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CRAFTEK  5 

SUMMARY  OF  RESULTS.  CONCLUSIONS 
AND  RECOMMENDATIONS 

The  results  presented  in  the  proceeding  chapters  demonstrate 
the  variations  in  control  policy  and  also  in  the  methods  for 
obtaining  the  same  when  the  tanks-ln-series  model  Is  substituted 
for  the  conventional  completely  mixed  tray.   The  mathematical 
models  used  were  adequate  to  represent  the  phenomenon  to  be 
investigated  without  having  to  go  into  the  complexities  of  the 
hydrodynamics  or  the  energy  balance  of  the  system.  Basically 
two  types  of  problems  have  been  treated,  the  time  optimal  problem, 
and  the  minimum  deviation  problem. 

In  Chapter  2,  the  reference  system  (conventional  completely 
mixed  single  tank  model)  was  set  up  and  the  time  optimal  problem 
was  treated.  The  control  policy  was  determined  so  as  to  transfer 
the  system  from  soue  initial  state  to  a  desired  final  state  in 
the  shortest  possible  time.  The  optimal  control  as  determined 
via  the  Maximum  Principle  was  found  to  be  purely  bang-bang  and 
the  results  are  presented  in  Table  1  and  Figures  15,  16,  1?  and 
18. 

In  Chapter  3«  the  two  tanks  in  series  system  was  set  up 
and  the  time  optimal  problen  i;as  once  again  treated.   Herein 
was  noticed  the  first  change  in  control  policy.   It  was  found 
that  the  final  desired  state  could  not  be  attained  by  merely 
using  bang-bang  control  and  hence  singular  control  was  resorted 
to.  The  results  are  presented  In  Table  2  and  Figure  31,  32,  33, 
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and  3^>   Also  as  a  check  on  the  berg-bang  policy  obtained  In 
Chapter  2  with  regard  to  the  reference  system,  the  time  optimal 
problem  of  Chapter  2  was  solved  again  by  using  singular  control. 
The  comparison  of  result;:  from  the  two  modes  of  control  is  made 
In  Table  3  and  Figures  36,  37,  38  and  39.  from  which  it  was 
quite  evident  that  for  the  reference  system,  the  bang-bang  policy 
indeed  achieves  time  optimality . 

Next  the  minimum  deviation  problem  was  considered  in 
Chapter  k   for  both  systems.   For  the  reference  system  it  seemed 
quite  obvious  that  bang-bang  control  should  suffice  once  more, 
as  it  did  far  the  time  optimal  problem,  however,  the  results 
turned  out  quite  opposite  to  the  obvious,  and  in  truth  singular 
control  had  to  be  used  to  obtain  minimum  deviation  optimality. 
This  was  confirmed  by  comparison  with  the  bang-bang  response  in 
Table  4. 

Perhaps  it  should  be  pointed  out  at  this  stage  that  a 
significant  observation  made  from  the  results  of  this  thesis  is 
the  different  conditions  under  which  the  necessity  for  singular 
control  should  be  recognized.   The  usual  conditions  is  that  the 
switching  function  becomes  zero  over. a  positive  interval  of  time, 
however,  as  pointed  out  by  Johnson  [12],  it  is  not  always 
possible  to  detect  this  when  ordinary  computer  searching  methods 
are  used.   Another  condition  is  that  sometimes  the  desired  final 
conditions  on  the  state  variables  cannot  be  attained  by  merely 
using  bang-bang  control.   When  such  a  situation  occurs  it  is 
only  natural  that  singular  control  be  used  to  achieve  the  final 
state  and  hence  perhaps  optimality.   This  was  the  situation 
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encountered  in  Chapter  3  vv  1  fch  the  two  tanks  in  series  system 
and  the  time  optimal  problem.  In  the  case  where  the  final 
state  can  be  arrived  at  by  bang-bcmg  control,  a  check  has  to  be 
made  on  the  Hamlltonian  to  see  whether  it  has  stayed  at  its 
minimum  throughout  the  transient  response  (In  the  case  where 
the  final  time  is  unspecified,  this  minimum  is  zero).   If  the 
Hamiltonlan  maintains  a  non  zero  value  then  this  is  an  indication 
that  even  though  the  bang-bang  policy  gets  us  to  where  we  want 
to  go,  in  doing  so  it  does  not  achieve  the  desired  optimality 
and  thus  singular  control  has  to  be  used.  This  precisely  is 
what  happened  In  Chapter  k   when  the  minimum  deviation  problem 
was  being  considered  with  the  reference  system.   Also  the 
superiority  of  the  singular  policy  over  the  bang-bang  policy  in 
achieving  a  minimum  deviation  is  shown  in  Table  4. 

In  the  latter  part  of  Chapter  k,   the  minimum  deviation 
problem  is  applied  to  the  mixing  pool  system  and  results  are 
tabulated  in  Table  5  and  corresponding  responses  are  shown  in 
Figures  k5,    k& ,   k7   and  W3 .  A  singular  policy  was  used  to  achieve 
a  minimum  deviation. 

In  Tables  6  and  7  are  given  the  comparisons  of  response 
times  of  the  reference  system  and  the  two  tanks  in  series  system, 
for  the  time  optimal  and  minimum  deviation  problems  respectively. 
Only  responses  from  singular  control  policies  are  considered. 
From  Table  6  we  see  that  fcr  the  tine  optimal  problem  the  refer- 
ence system  has  a  slower  optimal  response  time  than  the  two 
tanks  in  series  system  while  from  Table  7  we  see  that  for  the 
minimum  deviation  problem  the  reference  system  has  a  faster 
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TABLE  6 


COMPARISON  OF  RESPONSE  TIMES  TO  SINGULAR  CONTROL 
POLICIES  AS  APPLIED  TO  THE  REFERENCE  SYSTEM  AND  THE 
TWO  TANKS  111  SERIES  SYSTEM  (TIME  OPTIMAL  PROBLEM) 


Reference  System      Two  tanks  in  series 

system 
T,  mln  T,  mln 


Case  1  4.85  3.45 

Case  2  5.45  4.15 

Case  3  8.46  7,35 

Case  4  4.30  2.71 
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TABLE  ? 


COMPARISON  OP  RESPONSE  TIME  TO  SINGULAR  CONTROL 

POLICIES  AS  APPLIED  TO  THE  REFERENCE  SYSTEM  AND  THE 

TWO  TANKS  IN  SERIES  SYSTEM  (MINIMUM  DEVIATION  PROBLEM), 


Reference  System 

Two 

tanks  in  series 
system 

T,  min 

T,  min 

Case  1. 

7.70 

9-50 

Case  2 

7.80 

9.6c 

Case  3 

9-90 

10.95 

Case  h 

7.60 

9.^0 

' 
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optimal  response  time  than  the  two  tanks  in  series  system. 
These  results  from  Tables  6  and  7  show  how  the  response  times  of 
a  particular  system  can  be?  quickened  or  slowed  down  due  to  a 
change  in  the  objective  function,  i.e.,  dependence  of  response 
time  on  variation  in  objective  function.  Similarly  other 
dependences  may  be  noted.   For  Instance  in  the  time  optimal 
problem,  in  going  from  the  reference  system  to  the  two  tanks  in 
series  system,  the  control  policy  changed  from  bang-bang  to 
singular  thus  indicating  the  dependence  of  control  policy  on 
variation  in  model.   Also  in  the  reference  system,  when  going 
from  the  time  optimal  problem  to  the  minimum  deviation  problem, 
the  control  policy  changes  from  bang-bang  to  singular  which 
shows  the  dependence  of  the  control  policy  on  changes  in  objective 
function  when  the  model  is  held  constant.   The  problem  of 
meaningfully  e\'aluatlng  these  dependences  would  immediately  give 
rise  to  a  host  of  sensitivity  problems.   It  may  seem  vague  at 
this  stage  to  try  to  obtain  expressions  for  the  above  sensitivities 
by  the  generally  accepted  perturbation  techniques  that  have  been 
used  so  far  in  the  automatic  control  literature.   However,  as 
Srldhar  et  al.  [14]  have  put  it,  "Sensitivity  is  a  concept 
which  can  be  meaningfully  defined  only  by  considering  specific 
systems  and  their  particular  purposes  for  existence".   We  do  not 
wish  to  go  into  this  problem  presently  but  rather  point  out  the 
possibilities  for  further  research  in  this  area. 

To  some  readers  the  method  of  Grethlein  and  Lapidus  [15] 
may  not  seem  to  be  a  competent  method  to  arrive  at  a  singular 
policy  but,  as  was  stated  previously,  there  Is  no  general  method 
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at  present  by  which  a  singular  policy  can  be  decreed  optimal. 
Because  of  the  difficulty  in  obtaining  solutions,  the  present  use 
of  optimization  theory  for  control  systems  engineering  Is 
limited  to  determining  the  nature  of  the  optimal  control  law. 
This  information  is  useful  in  evaluating  alternate  control  schemes 
and  in  indicating  approximations  to  the  optimal  control  which 
may  be  more  easily  implemented.   On  a  similar  note  the  method  of 
Grethlein  and  Lapidus  was  incorporated  in  the  present  work  not- 
only  to  achieve  simplicity  in  obtaining  a  solution  but  also 
because  of  the  proven  versatility  of  the  method  in  a  somewhat 
similar  problem  as  described  in  [l5]« 

The  above  results,  discussion  and  observations  hold  only 
for  a  particular  model  with  a  particular  set  of  parameters, 
disturbance  and  initial  state.   It  should  be  recognized  that 
changes  in  the  various  parameters  or  in  the  assumptions  of 
Chapters  2  and  3  would  lead  to  diverse  situations  with  diverse 
results  and  making  any  generalizations  at  this  stage  may  not 
be  prudent. 

By  way  of  extensions  to  the  present  study  several 
recommendations  can  be  put  forth  for  further  study.  The  question 
of  sensitivity  analysis  has  already  been  raised  in  an  earlier 
paragraph.   In  addition  different  objective  functions  could  bs 
investigated  and  especially  of  the  type 

S  =   ;  [1  +  (x*  -  x(t))2]dt 
0 

wherein  both  time  and  deviation  are  minimised.   Also  various 
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modifications  can  be  made  on  the  mixing  pool  model.  In  Chapter 
3  the  assumption  was  made  that  the  vapor  rising  from  the  bottom 
tray  divided  into  half  and  each  half  entered  each  pool.   It 
would  be  Interesting  to  study  the  effect  of  uneven  division  of 
the  vapor  stream,  on  the  performance  of  the  model.   Another 
problem,  of   a  more  complicated  nature,  is  one  in  which  there  is 
feedback  between  the  mixing  pools. 

As  stated  above,  extensive  work  has  to  be  done  on  the  present 
model  before  any  general  conclusions  can  be  drawn.   However,  if 
the  present  Investigation  Induces  the  reader  to  pursue  the 
subject  further,  the  purpose  of  this  thesis  will  have  been 
fulfilled. 
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NOMENCLATURE 

B         Bottoms  product  flow  rate,  lb  mole/min. 

c  Constant  used  in  vapor-liquid  equilibrium  expression, 

y  =  mx  +  c. 

D  Overheads  product  flow  rate,  lb  mole/min. 

H  Hamlltonian  function. 

H  Bottoms  tray,  liquid  holdup,  lb  moles. 

B 
H  Top  tray  liquid  holdup,  lb  moles. 

L         Liquid  flow  rate  in  column,  lb  mcles/inin. 

L'         Reflux  flow  rate  which  maintains  desired  steady  state 
in  the  two  tanks  in  series  system,  lb  mole/min. 

m  Constant  used  in  vapor-liquid  equilibrium  expression, 

y  =  nx  +  C . 

S  ,  S       Initial  and  terminal  states  of  the  reference  system. 

1  2 

S' ,   S*      Initial  and  terminal  states  of  the  two  tanks  in 
-*■   ^      series  system.       , 

T  Terminal  time,  minutes. 

V  Vapor  flow  rate  in  column,  lb  moles/min. 

x  Light  component  composition  in  liquid  phase  on  top 

■*         tray  of  reference  system. Light  component  composition 
in  liquid  phase  in  pool  1  of  top  tray  in  two  tanks 
in  series  system. 

x  Light  component  composition  in  liquid  phase  on 

2  bottom  tray  of  reference  system.  Light  component 
composition  in  liquid  phase  in  pool  2  of  top  tray 
in  two  tanks  in  series  system. 

x  Light  component  composition  in  liquid  phase  on 

3  bottom  tray  of  two  tanks  In  series  system. 

y.         Light  component  composition  in  vapor  phase  in 


1 


equilibrium  with  liquid  phase  of  composition  x, . 


2.         Adjoint  variable  corresponding  to  state  variable 
x  .  • 
i 
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APPENDIX  A 

ANALOG  AND  DIGITAL  COMPUTER  PROCEDURES 

In  this  section  attention  will  toe  f ooussed  on  the  analog 
computer  circuit  diagrams  and  the  digital  computer  programs  of 
the  various  simulations  and  computations  used  in  Chapters  2,  3 
and  *K 

1.   ANALOG  COMPUTER  SIMULATION  OP  THE  REFERENCE  SYSTEM. 

The  system  equations  used  in  the  simulation  of  the  reference 
system  are  those  shown  in  equation  (11)  namely, 


dx- 

'T  "T   "        "t 


/Lm  -  Vm  -  Ln„         Vm        .    Lc  II  t\ 

w  -  ( — s: )xi +  c  2 +  r:  {a"1) 


and 


,dx2  ffl    .     fVm  +   F  t   L  -   V>  F*F        Vc 

dt    =   "    Hfi         (-  HB   ;        '  •-'    2  "  -HB-  *  HB 


x,.  --  mxL    -i-  c 


*v  -  mxi 

The  above  equations  were  magnitude  scaled  according  to  the 
scale  factors  shown  in  Table  A-l.   The  corresponding  magnitude 
scaled  equations  are 

T  T       T 

-  [i2]  -   •-  £■  [ij]  *    (Vm  t  F  t  I-  -  ?>[      ]   _  jg£  +  g   (A-2) 
E         '  B  B  B 

Od]  *  m[Xl]  +   c 


TABLE  A-l 
SCALE  FACTORS  FOR  REFERENCE  SYSTEM 


x  1.0  1.0  1.0 


X 


12? 


Problem 

Expected 
Maximum 

1 

Scale 
Factor 

Computer 

Variable 

Expected  Nax. 

Variable 

X  1.0  1.0  1.0       [x  ] 


'■*J 


1.0  1.0  1.0       [x  ] 

1)  «■  dj 
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The  time  scale  factor  p,  defined  as  the  ratio  of  computer  tine 
to  the  independent  problem  variable  t,  was  taken  as  unity;  thus 
one  second  on  the  analog  computer  represents  one  minute  of 
problem  time.  The  scaled  circuit  diagram  for  the  reference 
system  is  shown  in  Figure  A-l. 

2.   ANALOG  COMPUTER  SIMULATION  OF  THE  MIXING  POOL  MODEL  SYSTEM. 
The  system  equations  used  in  the  simulation  of  the  mixing 
pool  model  system  are  those  shown  in  eqiiatlon  (32)  namely , 

ffl   /Lm  -  Vm  -  2Lv      Lm      Vm      ,2Ln 
dt  -  K  H      '   1  r  Hm  x2  +  H_  *3    H_ 


dxo 


a?  - "  I:  *i  *  ^h*--^  *2  -  £  *3  (A-3) 


and 

x  =  0.5(mx  +  mx  +  2C ) . 
D        12 

The  above  equations  were  magnitude  scaled  according  to  the 
scale  factors  shown  in  Table  A-2.  The  corresponding  magnitude 
scaled  equations  are 

r*  ->    /La  -  Vra  -  2Lsr   -i  ,  Lit,  r   -i  ,  Vm  r-   l  ,  2Lc 
T  T         T         T 

-  ^  -  - 1  [xi]  +  ^f3^  -  e  t>3]      (-w+) 

T  X  T 

Fx„  -  Vc 


[x3]  ,  i  [x2]   -    (V*  t  r  t  L  -  V)^]  +  _  F_ 
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TABLE  A- 2 


SCALE  FACTORS  FOR  MIXING  POOL  MODEL  SYSTEM 


Problem     Expected     1 Scale     Computer 

Variable     Maximum      Expected  Max.    Factor    Variable 


X, 


1.0  1.0  1.0  [x,] 

S.          1.0  1.0  1.0  [x,.] 

2  L  Z- 

l.o  i.o  l.o  [>  ] 

*D       i-o  i.o  l.o  [xd] 


x 


13J 
[xD]  =  0.2?[X;l]  +  0,22  [x2]  +  0,56 

The  time  scale  factor  was  again  taken  as  unity.  The  scaled 
circuit  diagrams  for  the  nixing  pool  model  system  is  shown  in 
Figure  A-2. 

3.  DIGITAL  COMPUTER  SOLUTION  OF'TrlE  REFERENCE  SYSTEM. 

The  canonical  equations  of  the  reference  system  comprising 
of  equations  (11) ,  {Ik),    (16)  and  (17),  were  solved  along  With 
the  boundary  conditions  of  equation  (21a)  by  a  trial  and  error 
method  outlined  in  Figure  5«.  The  various  symbols  used  in  the 
program  are  explained  In  Table  A-3  and  a  sample  program  is  given 
in  Table  A-4. 

4.  DIGITAL  COMPUTER  SOLUTION  OF  THE  MIXING  POOL  MODEL  SYSTEM. 

Due  to  the  singular  nature  of  the  control  policy  for  the 
mi  -lng  pool  model  system,  the  sampled  data  method  of  Grethlein 
and  Lapidus  was  applied  to  the  system  equations  (32) >  The 
symbols  used  in  the  digital  computer  program  according  to 
Figure  J6   are  explained  in  Table  A-5  and  a  sample  program  is 
given  in  Table  A-6.   This  program  had  to  be  run  by  double 
precision  due  to  the  very  small  numbers  that  were  Involved  in 
the  computations. 
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TABLE  A-3.   NOMENCLATUBE  FOR  TRIAL  AND  ERROR 

SOLUTION  OP  MAXIMUM  PRINCIPLE  EQUATIONS. 


Program  Symbol 
A 
AL 

AUX 
DERY  (NDIH) 

FCT 

H 

NDIM 

OUTP 

PRMT  (1) 
PRMT  (2) 
PRMT  (3) 

PBMT  (k) 

PRMT  (5) 


X 
Y(I) 


Designation 

Switching  function 

Control  variable  L 

Auxiliary  storage  array 

The  differential  equations  that  are 
integrated  by  subroutine  RKGS 

Name  of  external  subroutine  which  defines 
the  differential  equations  DERY  (NDIH) 

Hamlltonian 

Number  of  differential  equations  to  be 
integrated  by  subroutine  RKGS 

Name  of  external  subroutine  for  handling 
output  from  subroutine  RKGS 

Initial  time  of  integration 

Final  time  of  integration 

Increment  of  independent  variable  t, 
over  which  integration  is  performed 

Tolerated  error 

Parameter  used  to  terminate  subroutine 
RKGS  at  any  desired  output  point 

Decimal  counter  used  to  monitor  the 
printing  of  results  at  various  points  in 
time.   Initially  Q  is  zero  and  after  inte- 
gration commences  Q  is  incremented  by 
0.01  each  time  results  are  printed 

Independent  vai'iable  t 

Dependent  variable  x 
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TABLE  A-4.   COMPUTER  PROGRAM  FOR  TRIAL  AMD  ERROR  SOLUTION 
OF  MAXIMUM  PBINCirXS  EQUATIONS 

EXTERNAL  FCTtOUTP 

DIMENSION  PRMT15) , Y I  3 1 , DERY I 3 » i AUX i 8, 3) 

COMMON  A.L.Q 
1  FORMATUH  ,6HIHLF  =i  13) 

PRMT(1)=0. 

PRMT<2)=4. 

PRMT(3)=.0l 

PRMTI4)=.00Q1 

AC* 1.3 33 

Q=0.0 

ND!M=5 

DERY(1)=NDIM 

DERY(l)  =  l.F.O/DERY(l) 

DO  /  I=2tNDIM 
7  DERYU)=DERYU) 

Y( 1)30.6300 

Y(2)=--.2767 

YI3)=C.O 

Y<5)=-.03 

Y(4)»-tl.+t.l413*AL-.12l3)*Y(5) )/( .207l*Al-.2071] 

CALL  RK.GS(PRMT,Y,DERY,NOIM,II!LF,rCT,OUTP,AUX> 

IF(IHLF.GT.10JWRITE(3,4)IHLF 

CONTINUE 
1?.  STOP 

ENO 
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TABLE  A-h.      (CON?. ) 


SUBROUTINE  FCT (X.Y.DtBVI 

DIMENSION  VII) ,DERY( 11 

COMMON  AL,Q 

DERY(i)  =  .44*AL*YC)-.5865*Y(I}-AL*YtlX„5865*Y(2H-.5ft*AL 

OERY(2)  =  .'i*AL*Y(  1!  ! .0986*Y( 2 1-.4*AL*Y 12 )-0. I486 

DERY(3)=l. 

DtRY(A)=-.A4*AL*Y('!!l.5865*Y(4)tAL<,Y(4)  «-.  'i*AL*Y  (  5  ) 

DERY<5)=-.5865*Yt<.!-.09S6*Y<5)«-.4*AL*Y(5) 

RETURN 

END 


SUBROUTINE  OUTPU,Y,DERY,  IHLF, NDIM.PRMT I 
DIMENSION  PRMrUI.YI  1),DERY(1) 
COMMCN  AL,0 
20  F0RMATI5X.F6.3) 

15  FOXHftfllX,*H  T  =,E11.4,6HX(1)  -  ,  E  I 1  .4,  6HX  f  2  )  =  ,  El  1  .4,  6!iX  (  3)  ■ 
ll.<i,6Hm)  =  ,E11.4,6H?.(2)  =,E11.4,4HL  =,F6.3) 

IF(X-O)  99,16,16 

16  WRITE!  3,  15  IX,  (Yd),  1-1,5),  AS. 
A=.4't*Y(l)*Yl'i)-Y(l)*Y(4)*.56*-Yl41*.4*Ytl)*Y<5)-.4*Y(2)*Yl5) 
IF1A)  17,17,18 

17  Al=1.333 

CO  TC  101  I 

18   AL=>833 
101  Q=Q+0.G1 

H=(.44*ftL*Ylll-.5865*Ytl)-AL*YllJ  f»5865*Y(  2)  «•  ,56*AL)*Y(  41 
1+(.4*AL*Y11M.0986*Y12)-.4*AL*Y(2)-.1'.S6)*Y(5)  +  1. 
WRITE  13,20)  H 
GO  TO  99 
9   PRMT1 51=1.0 
99  RETURN 
END 
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TABLE  A-4.      (CONT.) 


SUBROUTINE    RKGStPRM .YtDERY.NOiN.lHLFtFCT.OUTP.AUX) 
DIMENSION    VII)  ,DERY(i),AUXC8tl  5  ,  A  I  4  1  ,B  I  4  )  ,C(<> )  ,  PRMT  t  1 } 
DO    I    I=l,NOIM 

1  AUX(8,I)  =  .06666667>S'O--:RYm 
X=PRM1 (1) 
XEND=PRHT(2) 

H=PRMT(3) 

PRMT(5)«0. 

CALL  FCT (X,Y,DERY) 
C 
C.      ERROR  TFST 

1F(H*(XEND-XI)38,37,2 
C 
C      PREPARATIONS  FOR  RUNGE-KUTTA  METHOD 

2  AU>=.5 
A(2)=. 2928932 
A(3)=l. 707107 
A(4)=. 1666667 
B(l)=2. 
B(2I=1. 
B(3)±l. 

BC»)  =  2. 

CI1)=.5 

C(2>=. 2928932 

C(3)»l. 707107 

CU)  =  .5 
C 
C      PREPARATIONS  OF  FIRST  RUNGE-KUTTA  STEP 

00  3  I'UNOIN 

AUX(1,I)=Y(1) 

AUX(2,I)=DERY( I > 

AUX(3, I  1=0. 

3  AUX(6,t)=0. 
IRFC=0 

H=H  +  H 

IHLF=-1 

ISTEP=0 

IEND=0 
C 
C 
C      START  OF  A  RUNGE-KUTTA  STEP 

4  IF(  (X*H-XFND)*H)7,6,5 

5  H=XEND-X 

6  IEND»1 
C 

C      RECORDING  OF  INITIAL  VALUES  OF  THIS  STEP 

7  CALL  OUTP(X,Y,DERY, IREC, NOI M.PRMT ) 
IF(PRMT(5)KO,8,40 

8  irEST=0 

9  ISTEP=ISTEP<-1 
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C 

c 

C      START  OF  INNERMOST  RUNGE-KUTTA  LOOP 
J=l 

10  AJ=A(J) 
BJ=B(JI 
CJ=CiJ) 

DO  11  1=1,NDIM 
R1=H*DERY< I  ) 
R2=AJ*(R1-BJ*AUX(6,I)> 
Y(I)=Y(I J*R2 
R2=R2+R2+R2 

11  AUX(6,I)=AUX(6,t HR2-CJ*R1 
IF(J-4)12, 15,15 

12  J=JH 
IFIJ-3)13,14,13 

13  X=X+.5*H 

14  CALL  FCT(X,Y,OERY) 
GOTO  10 

C      END  OF  INNERMOST  RUNGE-KUTTA  LOOP 

C 

C 

C      TEST  OF  ACCURACY 

15  IF<ITEST)16,I6,20 
C 

C      IN  CASE  ITEST=0  THERE  IS  NO  POSSIBILITY  FOR  TESTING  OF  ACCURA 

16  DO  17  1=1,NDIM 

17  AUX(',,I)  =  Y{I) 
ITEST=1 
ISTEP=ISTEP+ISTEP-2 

18  IHLF=IHLF+1 
X  =  X-H 
H=.5*H 

DO  19  I=1,NDIM 

Y(I)=AUX(1,I) 

DERY(I)=AUX(2,II 

19  AUX(6,I)=AUX(3,  I ! 
GOTO  9 

c 

C      IN  CASE  ITEST=1  TESTING  OF  ACCURACY  IS  POSSIBLE 

20  IM0D=lSTEP/2 

IF ( I  STEP- IHOD- I  MOD  121,23,21 

21  CALL  FCT(X,Y,DERYI 
DO  21    1--1.NDIM 
AUX(5,I)=Y(I) 

22  AUX(7,I)=!)FRY(  I  ) 
GOTO  9 

C 

C      COMPUTATION  OF  TEST  VALUE  DELT 

23  OELT=C. 
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DO    24    I=i,NDIM 

24  DELT=0ELT+AUX(8,I  KA3S  (  AUXl«i.  !J-Y( II) 
ir(DELT-PRMII't)  >2C,2G,25 

C 

C  ERROR    IS    TOO   GREAT 

25  !F(IHLF-10)26,36,36 

26  DO  27  I=1,NDIH 

27  AUX(4,l)=AUX(5,I  ! 
ISTEP=ISTEP< ISTEP-4 
X=X-H 

IEND=0 

GOTO  18 
C 

C      RESULT  VALUES  ARE  GOOD 
2C  CALL  FCUX.Y.DERY) 

DO  29  I-1,NDIM 

AUX(1,1)  =  YU) 

AUX(2,I)*DERYU) 

AUX(3,I>=AUX(6,I) 

Y(I)=AUX(5,t) 

29  DERYf F)=AUX(7» I) 

CALL    OUTP(X-H,Y,DERY,IHLF,NDlM.PRMT) 
lF(PRMTt5l 140,30,40 

30  DO    31    I=1,NDIM 
Ytf )=AUXIl,I  ) 

31  0ERYU)  =  AUX(2,I) 
1REC=IHLF 
IF(IEND)32,32,39 

C 

C      INCREMENT  GETS  DOUBLED 
'   32  IHLF=IHLF-l 

lSTEP=ISTEP/2 

H=H+H 

IF<IHLF)4,33,33 

33  IM0D«ISTEP/2 
IFUSTEP-IM00-IM0D)4,34,4 

34  IF(DELT-.02*PRMTI4I)3S,35,4 

35  IHLF=IHLF-1 
ISTEP^ISTEP/2 
H=H+H 

GOTO  4 
C 
C 
C      RETURNS  TO  CALLING  PROGRAM 

36  IHLF=11 

CALL    FCTCX,Y,DERY> 
COTO    39 

37  IHLF=12 
GOTO    39 

38  IIILF=13 
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39  CALL  OUTPlX,Y,DERY,IHLF,NDIH,PRMT) 

40  RETURN 
ENO 
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Program  Symbol 
AL(J) 

AP 

ALEM 

DERY  (NDIM) 

DSV(J) 

FCT 

f 

NN 
NDIM 

OUTP 

PRMT  (1) 

PRMT  (2) 

PRMT  (3) 

PRMT  (4) 
PRMT  (5) 


Designation 

Array  containing  five  values  of  control 

variable  L  with  values  between  L    and 
t  max 

Hiin 

Control  variable  L  used  in  conjunction 
with  the  differential  equations 

Minimizing  decision  variable  L  from  among 
the  5  specified  in  AL(J) 

The  differential  equations  that  are  inte- 
grated by  subroutine  RKGS 

Deviations  corresponding  to  the  control 
variables  from  array  AL(J) 

External  subroutine  wherein  the  differential 
equations  DERI  (NDIM)  are  defined 

Number  of  current  sampling  period 

Total  number  of  sampling  periods 

Number  of  differential  equations  to  be 
integrated  by  subroutine  RKGS 

External  subroutine  used  for  purposes  of 
handling  output  from  subroutine  RKGS 

Initial  time  of  integration  of  current 
sampling  period 

Final  time  of  integration  of  current 
sampling  period 

Increment  of  independent  variable  t,  over 
which  integration  is  performed 

Tolerated  error 

Parameter  used  for  terminating  subroutine 
RKGS  at  any  desired  output  point 
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Q 

SUM 

SUMT 

X 


Designation 

Decimal  counter  used  to  ensure  the  com- 
pletion of  Integration  over  entire. 
sampling  period  before  evaluating  DEV(J) 

Value  of  deviation  computed  from 
objective  function  before  being  stored 
into  DEV(J) 

Summation  of  minimum  deviations  over  all 
sampling  periods  upto  and  including 
present  sampling  period 

Independent  variable  t 

Dependent  variables  x 


TABLE  A-6.   COMPUTER  PROGRAM  FOR  THE  SAMPLED  DATA  METHOD 
OF  GRETHLEIN  AND  LAPIDUS 

IMPLICIT  REAL*8< A-H,0-Z) 

EXTERNAL  FCT.OUTP 
C 
C 

DIMENSION  PRMT(6»,Y(4) ,DERY ; h) , AUX ( B, 4) ,YS(3) , AL 1 10  I , DEV ( 5 ) 

COMMCN  AP.N.YS, J,M,I I , ALEM, YA, YB t YC, YD, AL , SUM! 
C 
C      READ  NUMERICAL  DAI  A  ***************************************<. 

c 

NDIMM 
PRMT«1)=0.0 
PRMT(2I=0.05 
PRMT<3>=0.010 
P«MT(4J=0.1 
YS<1)=C. 71399 
YS(2)=0. 60986 
YS(3)=0. 26247 
ALU)*. 833 
AL(2>=.953 
AL<3)=1.083 
ALI1I=1.208 
AL(5)=1.333 
N=l 

NN=25C 
11  =  1 
M-0 

SUMT=0„0 
C  / 

c 

C      READ  PARAMETERS  FOR  ITERATION  KITH  ALII)  ******************* 
C 

100  J=l 
AP=AL(J> 

IF(M.EQ.l)  GO  TO  150 
V«l)»C. 71467 
Y(2)=0. 60909 
Y(3)=0. 24899 
V(4)=0.0 

GO  TC  101 
150  Y( 11-YA 
Y(2)=YB 
YI3)=YC 
Y(4)=YD 

101  CALL    RK.GS(PRMT,Y,D£RY,NDlMttHLF,FCT,OUTP,AUXI 
CONTINUE 

C 

C      READ  PARAMETERS  FOR  ITERATION1  WITH  AL(?)  *******************> 

C 

AP=AL(J) 

IF(H.EQ.l)  GO  TO  160 
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Y(1)=0. 71467 

Y12>=0. 60909 
Y(3!=C. 24899 
YCVI'0.0 

GO  TC  102 
160  Y( l)  =  YA 
Y(2!=YB 
Y(3)=YC 
Y(4)=YD 

102  CALL  RKGS(PRMT,Y,DERY,NDIM, IHLF, FCT, OUTP , AUX ) 
CONTINUE 

READ  PARAMETERS  FOR  ITERATION  WITH  Al ( 3 )  ******************* 

AP=AL( J) 

IF(H.EQ.l)  GO  TO  170 
Y( 11=0,7146  7 
Y(2)=C. 60909 
Y(3>=0. 24899 
Y(4)=0.0 
GO  TC  103 
170  Y(l)~YA 
Y(2)=YR 
Y(3)=YC 
Y(4)=Y0 

103  CALL  RKGS(PRMT,Y,DERY,N01M, IHLF, FCT, OUTP, AUX) 
CONTINUE 

READ  PARAMETERS  FOR  ITERATION  WITH  AL(4)  ******************* 

AP=AL(J) 

IFIM.EQ.1)  GO  10  180 
Yl  U=C. 71467 
Y(2>=0. 60909 
Y(3)=0. 24899 
Y(4)  =  C0 
GO  TC  104 
ISO  Y( 1)=YA 
Y(?I=YB 
Y(3)=YC 
Y(4)=Y0 

104  CALL  RKGS!PRKT,Y,DFRY,NDIM, IHLF , FC I , OUTP, AUX I 
CONTINUE 

READ  PARAMETERS  FOR  ITERATION  WITH  AL ( 5 )  ******************* 

AP=AL(J! 

IF(K.EQ.l)  GO  TO  190 

Yd  )=C. 71467 

Y( 2) =0.60909 
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Y>3)=C. 24899 
Y(4)  =  C0 
GO  TC  105 
190  Y(1)=YA 
Y(21=YB 
Y13I=YC 
Y(A)=YD 

105  CALL    RKGSIPRr-IT.Y.DERYtNDIM.IHLF.FCT.OUTP.AUX) 
CONTINUE 

READ  PARAMETERS  FOR  ITERATION  WITH  ALEH   ******************* 

AP=ALEK 

IF(H.EO.l)  GO  TO  2C0 

Yll)»0.7K67 

Y<2)=0. 60909 
Y13)=0. 24899 
Y(4)=0.0 
GO  TC  106 
200  Y(1)=YA 
Y(2)=YB 
Y(3)-YC 
Y(4)=YD 

106  CALL  RKGSIPRMT,Y,0ERY,NDIK,IHLF,FC1,0UTP,AUX) 
CONTINUE 

RETURN  FOR  REITERATION  AT  (N+11TH  SAMPLING  PERIOD*********** 

II^l 

M=l 

N=N 

AN»N 

Y(l) 

Y(2) 

Y(3) 

VHI 

PRM1 

PRMT 

IF(N 

STOP 

END 


1 
-I 

=  YA 

=  YB 

=  YC 

=  YD 

<1)=AN*0.05 

(2l=PRHT(l)+0.05 

.LE.NNI  GO  TO  100 
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SUBROUTINE  RKGSIPRM  , Y , DERY , NO  I H , I FLF , F CT, CUTP, AIM  ) 
IMPLICIT  REAL*8IA-H,G-Z) 


01  HEN  SIGN  Y(1),DERY(  1)  , AUX ( 0, 1 ) , A < 4 ) ,8 < 4 )  ,  C (4  l.PRMTtl) 

X=PRPT(1> 

H  =  PRI>'T(3) 

PRMT(5)=0. 

CALL  FCT(X,Y,DERY> 


PREPARATION'S  FOR  RUNCE-KUT FA  METHOD 

2  A(l)=,5 
A(2)=. 2928932 
A(3)=l. 707107 
A(4)=. 1666667 
B(l)=2. 
B[2)=l. 
B(3)«l. 
B(4)--2. 
C!l)=.5 
C(2)=. 2928932 
C(3)=l. 707107 
C  ('.)  =  .  5 

PREPARATIONS  OF  FIRST  RUNGE-KUTTA  STEP 

00  3  t=l,NDIM 

AUXU.I  )  =  YiI)  / 

AUXI2.I )=DERYl  I  I 

AUX(3,I 1*0. 

3  AUX I  6, I 1-0, 

RECORDING  OF  INITIAL  VALUFS  OF  THIS  STEP 
7  CALL  CUTP(X,Y,DERY, I REC , NDI M, PRMT I 
IF(PRCT(5) 140,8,40 


START  OF  INNERMOST  RUNGE-KUTTA  LOOP 
8  J=l 

10  AJ=A(JI 
BJ=B(J) 
CJ=C(J) 

DO    11     I-l.NOIM 
Rl=H  +  [JERY( I J 
R2-AJ*(FU-FiJ*AUX(6,I  )  I 
Y(I)=Y( I I+K2 
R2=R2+R2+R2 

11  AUX(6,D=AUXI6,I  )*R2-CJ*R1 
IHJ-4112, 15,15 

12  J=J+1 
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irU-3J13,U,13 

13 

X=X+.5*H 

14 

CALL  FCT(X,Y,DER 
GOTC  10 

c 

END  CF  INNERMOST 

c 

c 

15 

DO  2$  [=l,NOIM 
AUX(I,I  )=YU) 

RUNGE--KUHA    LGQP 


AUX(2,I  i=t)ERY(  1  ) 

29  AUXlfci D-AUXI3, I ! 

CALL  CUTP!X,Y,DERY, IHLF.NDIM, PRMT) 
IFIPRMT(5)K0,30,40 

30  00  31  I=l,NOIM 
Y(  t  )  =  AUX(1, I) 

31  DERYI l)  =  AUX(2,  i! 
GO  TO  8 

'.0  RETURN 
END 


SUBROUTINE  FCT I  X , Y , DERY ) 
IMPLICI1  REAL*8(A-H,G-Z) 


DIMENSION  YUI  ,  DERYU)  , YS < 3) , AL I lu > 

COHMCN  AP,N,YS, J,M, II , ALEM . YA , YB , YC, YD, AL, SUMT 

DEFINE  DIFFERENTIAL  EQUATIONS   ****************************** 

DERYtl)=-l.'36*AP*Y{l)-.56652*Y(  I )  <• ,  ','t*AP*Y  (  2  )<• .  5B65  2*Y  (  3  1  H  .  1 
DERY (2 )^2.0*AP*Y(l)-.5e652*Y(2)-2.C*AP*Y(2)+.58652*Y(3) 
DERY(3»  =  AP*.4*Y(2)+O.Oy86*Y(3)-0.'.*AP*Y(3)-0.1't8'592 
DERY  (4  )* '( .  B5123-.22*Y  ( 1 )  -.22*Y  ( 2 )-.  56 )  * 
M.85123-.22*YU)-.22*Y<2)-.56) 
RETURN 
END 
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SUBRCU1 1NE  OUTP(X,v,DERY, IHLF,NDIM,PRMTJ 
IHPLICiT  REAL*S(A-M,0-Z) 
C 

c 

01  KENS  I  ON  PRMTie'sYCt)  tOERY(4),VS(3),AL(lO»,DEV(rj) 
COMMCN  AP,N,YS, J,M, I  I  ,  ALEM,  YA,  YB  ,  YC,  Y[),  AL  t  SUMT 

1  FORMAT  (5X,<iHCEV(  ,  II, Ml)  =  ,!-"].'.. B) 

2  FORMAT! 5X.4HN  =  ,  I3i 

3  FORMA T<5Xt5HAP  ■  ,F10.5) 

'.  FORMATI5X.7HXU)  =  ,  F10.  5,  5X,  7HX  (  2  I  =  ,  F10.5.5X  ,7HX(  3 )  -  ,F1C 

15X,7HX(4)  =  ,F14.9) 
5  FORMAT(5X,7H0EVM  -  ,F14.8I 
■111  BN=N 

Q»BN*0.05 
IF(X-C)  20,10,10 
20   RETURN 
10   IF(II.LE.O)  GO  TO  50 

r. 

C     CALCULATE  DEVIATIONS  FROM  STEADY  STATE  FOR  JTH  REFLUX  RATE****: 
C 

SUM=SUMT* 
1  (.28247-YI3) I *( .20247-YI 3) ) 

DEVU)»SUN 

WRITE(3,l)  J,  DEV(J) 

IFIJ.GE.5)  GO  TO  30 

J=J  +  1 

PRMT(5)=1.0 

RETURN  i 

C 

C      FIND  MINIMUM  DFVIATION   *****♦*********************«♦********' 
£■ 

30   0MIN=DEV(1) 

JMJN=J 

DO  AC  J=2,5 

IFIDEVI Jj.LT.DMIN)  JMIN=J 

IFIDEVI Jl.LT.DMIN)  DMIN=DEV(J) 
40   CONTINUE 

AL1=ALII) 

AL2=AL(2! 

AL3=AL(3! 

AL4=AL(4! 

AL5=AL(5I 

OEV]=CEV( 1) 

DEV2=CEV(2) 

DEV3=CEV(3] 

DEV4=0EVC41 

DEV5=DEV(0) 

IF! JMIN.EQ.l)  ALEM-ALI JMIMI 

IFUMIN.EO.l)  WRITE ( 3. 5) DEVI JM IN  I 

IFUMIN.EQ.l)  DEVJ=DEV(JMIN) 
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IFtJCIN.EQ.5)     ALEM=AL(JMIN) 

IFUHN.EQ.5I    WRITE<3f5)DEV(JMINI 

IFUPIN.EQ.5)    DEVJ=DEV(JMIN) 

IFt JMIN.EQ.2)  CALL  POL  1 < DE VI , PEV2 , CEV3, AL I , AL2» AL3, ALEM, DEV J 

IF(  JPIN.EQ.3)  CALL  POL  1  (  OEV^,  UL1  V3  ,  DEV4  t  AL2  ,  AL3  ,  AL't,  ALEM,  DEV  J 

IFIJPIN.E0.4)  CALL  P0Ll<DEV.3,r>EV4.CEV5,  AI.3,  AL4|AL5,ALEM,DEVJ 

AP«ALEM 

SUMT=DEVJ 

11=0 

PRMT(5)=1.0 

RETURN 
C 
C      WRITE  SOLUTIONS  FOR  NTH  SAMPLING  PERIOD******************** 


i 


50   WRITF;(3,2)  N 
WRITE(3,3)  AP 

WRITE13.41  Y(1),Y(2),Y(3),Y(4) 
YA=Y(1) 
YB=YI2) 
YC=Y(3) 
YD=Y('i> 
PRHT(5)=1.0 
RETURN 
ENO 


SUBROUTINE  POL l( Fl , F2t F3, Rl ,R2, R3.RMIN, FMIN) 

IMPLICIT  RFAL*R(A-H,0-Z) 
200  F0RPATC>X,7HRMIN  =  ,F10.5f 5X, 7HDEVM  =  .Fl't.S) 
C 

C      CALCULATE  COEFFCIENTS  OF  3  POINT  LAGRANGI AN***************** 
C 

112  DI=CR1-R2)*IR1-R3) 

D2=(R2-R1 )*IR2-R3) 

03=(R3-R1)*[R3-R2) 

EI=F1/D1 
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E2=F2/02 

E3=F3/03 

GUR2  +  R3 

G2=RHR3 

G3=R1+R2 

H1=R2*R3 

H2=Rl*R3 

H3=R1*R2 

A=E1+E2«E3 

B=E1*G1+E2*G2+E3*G3 

C=E1*H1+E2*H2+E3*H3 
C 

C  CALCULATE    RMIN    AND    FMIN 

C 

RMIN=B/(2.0*AI 

FKIN=A*RMIN*RMIN-B*RMIN+C 

WRITE(3,20C)  RMIN, FMIN 

RETURN 

END 
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APPENDIX  B 

STEADY  STATE  AND  LIMITING  CASE  ANALYSES 

In  order  to  determine  the  physically  realizable  bounds  of 
the  two  systems  studied  in  this  work,  steady  state  and.  limiting 
case  analyses  have  been  carried  out  for  each  system.  This 
Appendix  deals  with  these  analyses  starting  with  the  reference 
system  analysis. 

1.   REFERENCE  SYSTEM  ANALYSIS. 

The  performance  equations  of  the  reference  system  are 


dx 

m 


5-.   (La.-  Via  -  L)      Vm 

*1   +  U    Xo  +  TJ~ 


Lc 

H  ^'H 
T       T 


dx 


Px. 


— 2  _  _L  x  _  .(Vm  H-  F  +  L  -  V)„   ,  »] 
dt  -  H  *1  H        X2  +  T~  ~   H~ 

B  B     '        B     B 


(B-l) 


(B-2) 


Multiplying  equations  (B-l)  and  (B-2)  by  H  and  H  respectively, 

T       B 

and  setting  the  derivatives  with  respect  to  time,  of  the 
resulting  equations,  equal  to  zero  yields 


(Lm  -  Vm  -•  L)x  +  Vm  x  =  -  Lc 
L  x.  -  (Vra  +  B)x  =  Vc  -  Fx 


which  can  be  put  in  matrix  form  as 

[Lm  -  Vm  -  L     Vra 
L      -(Vm  +  B). 
Equation  (B-5)  is  of  the  form 


r 


l.  X 


-Lc 

Vc  -  Fx 


(B-3) 
(B-'t) 


(B-5) 
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AX  =  B 


(B-6) 


from  which  X  can  be  solved  for  uniquely  by  applying  Cramer's 
Rule,  provided  the  coefficient  matrix  A  is  nonslngular.  Hence, 
the  system  has  the  unique  solution 


-Lc 

Vm 

Vc  -   Fx„ 
F 

-(Vm  +   B) 

xl  - 

Lm  -  Vra  - 

L 

Vm 

1       * 

■•(Vm  +   B) 

(B-V) 


and 


X2  = 


Lm  -  Vm  -  L 


-Lc 


Vc  -  Fx 


Lm  -  Vm  -  L 


Vm 


L        -(Vm  -i-  B) 

Equations  (B-7)  and  (B-8)  can  be  further  simplified  to 

2 

LVmc  +   BLc   -   V  'mc  +   FVmx,? 

V  m     -   LVm     +   BVra   -   BI,m  +   BL 


and 


LVcc  -  V2mc  -  LVc  -  FLex„  +  FVmxT,  +  FLx„  +  L2 
,  F F      F 

V2m2  -  LVm2  +  BVm  -  BLm  +  BL 


(b-8) 


(B-9) 


(B-10) 


Hence  x_  and  x  can  be  generally  written  as 


152 


x  =  x  (L,  V,  P,  B,  m,  o,  x?) 

(B-ll) 
xg  =  x2(L,  V,  F,  B,  m,    o,  x  ) 

however,  in  find  o  are  fixed  by  the  equilibrium  relation  of 
equation  (4a) 

y,  =  0.44  x.  4  0.56     1  =  1,  2 

and  B  is  itself  a  function  F,  V  and  L  as  can  be  seen  in  equation 
(la), 

B  =  F  +  L  ••  V 

whereby 

xl  =  xi^L'  Vl  F<  XF' 

(B-12) 
x2  -  x2(L,  V,  F,  xF) 

Thus  x  and  x  have  been  reduced  to  being  dependent  on  the  four 

parameters  L,  V,  F  and  x  .   It  is  obvious  that  for  each  set  of 

F 
values  for  the  above  parameters,  x  and  x  have  unique  values 

corresponding  to  a  unique  steady  state;  and  also  innumerable 
combinations  of  the  above  parameters  will,  give  rise  to  innumer- 
able steady  states.  Hence,  in  order  to  narrow  down  the  range 
of  steady  state  solutions  of  the  reference  system,  we  shall 

restrict  ourselves  to  the  case  where  V,  F,  and  x  are  held 

F 

constant  as  prescribed  by  equation  (4),  i.e. 

V  =  1.333  lb.  moles/min. 

F  r:  0.5  lb.  moles/min,  (B-13) 
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x  =  0.65  (B-13  cunt.) 

F 

Utilizing  the  overall  material  balance  and  a  material  balance 

around  the  reflux  drum,  >:e  have 

F  =,  D  +  B 
and 

V=L  +  D 

respectively;  and  using  the  values  of  V  and  F  as  in  equation 

(B-13)  we  find  that  L  is  restricted  to  values  between  0.8 33  and 

1.333«  Physically  this  means  that  the  overhead  reflus  flow 

rate  can  be  varied  from  a  minimum  of  O.833  lb.  moles/ruin,  to  a 

maximum  of  1.333  lb.  moles/min. 

The  steady  state  analysis  consists  of  using  values  of  V,  F 

and  x  given  by  equation  (B-13),  in  equations  (B-9)  and  (B-10) 

and  evaluating  x  and  x  as  L  is  varied  between  its  lower  and 
1      2 

upper  limits.  This  same  procedure  is  then  repeated  for  different 

values  of  x  (since  ultimately  we  wish  to  determine  an  optimal 

r 

path  for  L  so  as  to  counteract  a  known  disturbance  in  x„). 

r 

The  results  are  plotted  in  figure  B-l,  for  x  =  0.55,  O.65, 

F 

and  0.75.  Also  in  the  plot  are  Included  curves  for  the  overhead 
composition  x  since 

x  =  jr.  «=  mx^   +  c 

Figure  B-l  acts  as  an  a:5d  in  determining  the  range  within 
which  the  reference  system  can  operate  for  the  given  set  of 
parameters. 
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Fig.B-l.      L  v/s   x    with    Xf  cs     parameter     for 
reference     system. 


2.   MIXING  POOL  MODEL  SYSTEM  ANALYSIS 

The  pcrf ormance  equations  of  the  mixing  pool  model  ara, 


dx 

"dt 


1   (Lm  ~  Vm  -  ?L)  _   .  Lm  _   .Vm  „   .  2Lo 


"•1  T  H 


-2  -r  g   -3 

T  T 


dx 


2  _  2L    _  (Vm  +  2L) 


dt  -  H   A 


Va 

x2  +  ^  X;j 


Fxr 


5^2  _  -L  x  -  1Y.MJ:  F  +  L  -jrj  "J  _   Vc 

dt  -  HB  X2  Hfi         3  '  HB    HB 


(B-HO 
(B-15) 
(B-16) 


Multiplying  equations  (B-1'0  and  (B-IS)  by  H  and  equation 

T 

(B-l6)  by  Hn  and  setting  the  time  derivatives  in  the  resulting 
equations  equal  to  zero  gives 


(Lm  -  2L  -  Vm)x,  +  Lmx?  +  Vmx^  a  -2Lo 


2Lx1  «  (Vm  4  2L)x2  +  Vmx-,  ~   0 


Lx„  -  (Vm  +  F  +  L  -  V)x„  a   Vo  -  Fx,, 


(B--X7) 
(B-18) 
(B-19) 


which  in  matrix  form  becomes 


Lm  -  2L  -  Vm      Lm  Vm 

2L      -(Vm  +  2L)  Vm 

0  -(Vm  +  F  +  L  -  V) 

-  2Lc  " 
0 
iVc  -  Ft 


~3J 


(B-20) 


Again,  equation  (B-20)  is  of  the  form 
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AX  =  B 


froiii  which  X  can  be  solved  for  by  using  Cramer's  Rule,  provided 
the  coefficient  matrix  A  is  nonslngular.   Hence,  the  unique 
solution  of  the  system  is 

Vm 
Vm 
•-(Vm+F+L-V) 


-2Lc 

Lm 

0 

-(2L+Vm) 

Vc-Fx„ 
F 

L 

SI  - 

A 

Lm~Vm-2L 

-2Lc 

2L 

0 

0 

Vc-Fx^ 

Vm 

Vm 

-(Vm+F+L-V) 


(B-21) 


(B-22) 


and 


' 

Lm-Vm~2I. 

Lm 

-2Lc 

2L 

-(2L+Vm) 

0 

0 

L 

Vc-Fx 
F 

3  _ 

(B-23) 


where 


Lm-2L-Vm 

Lm 

2L 

~(Vm+2L) 

0 

L 

V'"i 

Vm 
•(Vm+F+L-V) 


15? 


Equations  (B-21),  (B-22),  and  (B--2.")}  can  be  further  simplified 

to 

K 

^  -  ^  (3-2/1) 

where 


N,  r-  -4cL  -  kFol,     +  (')-c  -  4mc)L  V  -  (2Fmc  +  ft  i  +  2Fnx  ) 
1  P       F 


where 


LV  +  (W.  -  m2o)LV2  -  Fm2x  V2  +  m2cV3 


D  =  ('to  -  'J) I,3  +  (tea  -  4F)L2  +  ('to2  -  8m  +  4)L2V 

+    (Pm2  -  tem)LV  +    ('Ma  -  5m2  +  m3)LV2  -  Fin2V2  -   (m3-  ro3)V3 

N? 
x2  =  -f  (B-25) 


N     b   --4el.3  -   teiL2  +    (4o   -'  *tonc)L2V  +    (FxFm2  -  ^FZjajLV 
+    (4-mc  -  n^cjLV^  -   Fxja2V2  *   m2cV3 


N 

*3  =  -g  (B-26) 

where 

N     =   -if-cL3  +    (tex„m  -  texjL2  +    (^0  -   4ffie)I,2V 

+   (PXjja  -  tetptajLV  +    (tone  -  m2c)LV2-  Fm2xFV2+  m2oV3 

Hence,   x   ,   x     and  x     can  be  generally  written  as 
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Xj  =  x1(L,  ,V,   F,   m,   c,   af) 

x^  =  x2(L,   V,   F,   m,   c,   xp)  (B-27) 

x     a  x  (L,    V,    F,   m,    c,    x   ) 

which  can  be  further  reduced  as  in  the  reference  system  to 

x    =  x_(L,   V,   F,   xj 
11  F 

X     =   x   (L,    V,    F,    x   )  (B-28) 

*Z  ti  r 

x     e   x   (L,    V,    F,    x_) 
3        3  F 

Arguing  along  the  same  lines  as  in  the  reference  system  analysis, 

we  have  V,  F,  and  x„  held  constant  as  in  equation  (B-13)  and 

r 

find  that  L  is  restricted  to  taking  on  values  between  0.833  and 
1.333-   Hence,  the  steady  state  analysis  consists  of  holding 

V,  F,  and  x  constant  in  equations  (B-2^),  (B-25)  and  (B-26) 

F 

and  evaluating  x,  ,  x  ,  and  x  as  L  is  varied  between  its  lower 
and  upper  lomltsj  the  procedure  being  repeated  for  different 

values  of  x  • 

F 

The  results  are  plotted  in  figure  B-2  for  x  =  0.55,  0.65, 

r 

and  0.?5«  Also  in  the  plot  are  included  curves  for  the  overhead 

composition  x  since 
D 

x  =  0.22(x,  +  x  )  +  C 
D        12 

Figure  B-2  acts  as  an  aid  in  determining  the  range  within 
Which  the  mixing  pool  model  can  operate  for  the  given  set  of 
parameters. 

Diagrams  of  the  type  of  figures  B-l  and  B-2  should  form  an 
integral  part  of  any  sort  cf  analysis  work  involving  systems 
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Fig.B-2.     L  v/s  x  for   tanks-in-series    model    using 
xf  as    parameter. 
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with  several  parameters  since  they  give  at  a  glance  an  overall 
picture  of  the  limits  within  which  the  system  operates. 


l6l 


append;;-';,  c 


COMPUTATIONAL  ALGORITHM  OF  POHTHXAGIN ' S  MAXIMUM  PRINCIPLE 
AND  STATEMENTS  OF  THEOBEMS  ON  EXISTENCE  OP  OPTIMAL  CONTROLS 


A  continuous  process  is  best  represented  by  a  set  of 
simultaneous  1'irst-order  differential  equations;  the  number  of 
equations  depending  on  the  order  of  the  system.   In  general, 
an  nth  order  continuous  system  can  be  represented  by  n 
simultaneous  first  order  differential  equations  by  using  the 
concept  of  state  space.  The  differential  equations  are  called 
the  performance  equations  of  the  process. 

We  shall  now  proceed  to  state  the  algorithm  based  en  the 
Maximum  Principle  associated  with  continuous,  autonomous  (i.e. 
time  does  not  appear  explicitly)  systems. 

C-l.   STATEMENT  OF  THE  MAXIMUM  PRINCIPLE  ALGORITHM. 

Let  the  performance  equations  of  a  process  have  the  form 

dx 

-gi=,  f^x-Jt),  x2(t),  ....  xn(t),  e^t),  c2(t),  .,.,  ep(t)) 


(C-l) 


or  the  vector  form 


ff  =  f(x(t);  e(t))  (C-2) 


where  x(t)  is  an  n-dimenslonal  vector  function  representing  the 
state  of  the  process  at  time  t  and  6(t)  is  an  r-dimensional 
vector  f unction  representing  the  decision  at  time  t, 

(In  the  case  where  the  system  is  linear,  the  equation 


162 


comes 

dXl 

dt    * 

n 

r 

-     £ 

j=l 

aij  xj  +  ^  bik  ! 

>k) 

(C-3) 


with  Initial  conditions 


x(t)     =  x  (0)  =  a       1  =  1,2 n        (C-k) 

1    t=0    1       1 


The  problem  associated  with  such  a  process  is  to  find  or  choose 
0(t)  subject  to  the  constraint 

6(t).£§  (C-5) 

(where  (p  is  an  r~dimensional  space),  which  makes  a  function  of 

the  final  values  of  the  state 

n 
S  =  X  c,  x,  (T),       C,  =  const.  (C-6) 

1=1  i  x  1 

a  maximum  (or  a  minimum).  The  function  S  to  be  extroniized  is 

called  the  objective  function,  and  the  decision  vector  function 

so  chosen  is  called  the  optimal  decision  0(t). 

The  first  step  in  solving  such  a  problem  is  to  Introduce  a 

new  state  variable  so  that  the  desired  form  of  the  objective 

function  may "be  obtained.   For  example,  if 

T 
/ 
0 

is  the  objective  function,  then  an  additional  state  variable 


S  =   /  (1  t  82)dt 


x    is  introduced  such  that 
n+3 
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T 
x   (T)  =  S  „  /  (1  +  e")dt 


-&~m  X  *   e  •      xml(,)  ■  ° 

and  hence  the  objective  function  novi  becomes 

S  .  X    (T)  (C-?) 

n+l 

Finally  an  n-dimenslonal  adjoint  vector  f:(t)  Is  Introduced  along 

with  a  Hamiltonlan  function  H  which  satisfy  the  following 

relations, 

n 
H(z(t),  x(t),  6(t)  =  2  z.    f(x(t),  6(t)  (C-8) 

1.1  x     * 


and 


where 


£fi=.iiL=_  I  z.ih 


--*$   "jflzjpr     I- !.*.•••.■    <c-9> 


The  set  of  equations  (C-l)  and  (C-9)  constitutes  a  split 
boundary  value  problem,  whose  solution  depends  on  6(t).  The 
optimal  decision  vector  function  ¥(t)  which  makes  S  an  extremua 
also  makes  the  Hamiltonlan  an  extremum  for  all  t,  0  <  t  <  T. 
Furthermore,  the  maximum  (or  minimum)  value  of  H  is  a  constant 
for  every  t. 

C-2.   DERIVATION  OF  THE  ALGORITHM 

In  this  section  we  shall  consider  3  basic  types  of  probleiss 
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and  derive  the  Maximum  riinciple  algorithm  In  each  case.  For 
the  case  where  the  final  time  is  fixed  we  have  two  problems,  a 
fixed  right-hand  and  a  free  right-hand  problem,  depending  on 
whether  the  final  conditions  on  the  state  variables  are  given  or 
not,   The  third  problem  arises  for  the  case  where  the  final  time 
Is  left  unfixed. 

C-2.1.   PINAL  TIME  FIXED  AND  FREE  EIGHT-HAND  PROBLEM 

The  system  will  be  the  same  as  governed  by  equations  (C-l) 
and  (C-4).   The  objective  function  is  defined  as  in  equation 
(C-6),  and  is  to  be  maximized.  Next  the  adjoint  vector  function 
along  with  the  Hamiltonian  is  introduced  as  defined  In  equations 
(C-8)  and  (C-9). 

Now,  if  6(t)  represents  the  optimal  decision  vector  such 
that  S  attains  its  maximum  value,  then  a  small  perturbation 
6  0(t)  will  move  the  optimal  state  vector  x(t)  by  a  small 
deviation  5x(t),  (see  figure  C-l).   That  is,  if 

dx. 

-^-  =  f^xU);  e(t))  (CIO) 

represents  the  optimal  system,  then  after  perturbation  we  have 


££  (*j  +  6Xl)  =  f  (x  +  6x;  8  +  60)  (Oil) 

Upon  subtracting  equation  (C~10)  from  equation  (C-ll)  we  have 

dT  6xi  "-:  fl(*  *  6x'  «  H-  60)  ~  ft(xt   ~)  (C-12) 


Multiplying  both  sides  by  the  sdjoint  variables  z   (t)  and  then 


,M 


t=0 


Fig.  C- 


■ 


-».  ■«,  - 


Perturbation     trajectory    due    to    change    in 
decision     variable    for    two    dimensional    sys 
tern.     Notice     how     £x(O)  =  0     when     initial 
values     are     preassigned. 
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taking  the  sum  over  the  subscript  1  yields , 

n  n 

Z    z  7%t>x>    =  £  z,  f,(x  +  6xs  5"  +  6B  )  -  f,(x|  6  )      (C-13) 
1=1  iat  1   1=1 


Next  Integrating  both  sides  of  equation  (C-13)  with  respect  to 
time  from  t  =  0  to  t  =  T  gives 


T  n 


d 


/   £  Z,i5r.dt  =  /   £  z   f  (x  +  6x;  6  +  60  )  -  f.  (xi  0  )at 


0   1=1 


'IcfT'-i 


0   1=1 


1  i 


(C-UO 


Now  consider  the  differentiation  of  £  z,  6x 


1=1 


1        1 


d      n  n  dx,  n  . 

dt  Azi  6xi  -  /  -dT  6xi  *  .^zi  It  6xi 


i=l 


i=l 


1=1 


(c-15) 


hence 


T  n    j,         n 

/   £  z     yrr   6x  dt  =  £  z   t>x 
0  1=1  *  dt   1     1=1  3   l 


T     T  n  dz 

-  /  £  -1-  6x-  dt     (C-16) 
0   0  1=1  ^ 


dt    1 


Therefore  combining  equations  (C-l^)  and  (C-16)  we  get 


£  z,6x 
i=l 


T  n  dz  T  n 

/  z   -~6x  dt  t   /  £  z   f  (x  +  6xi  8  +  66) 


1   *  lo  "  0   1=1  "   1    ''  0J    1=a  i   iV 


fji(x(  0)dt 


(C-17) 


Since  the  initial  state  values  have  been  preassigned  in  equation 
(C-4) 


6x1(0)  =  0 


(see  figure  C-l) 


and  also  since  z  (T)  =  C,  ,  we  have 
i       1 


n 

zi   6xi 


n=l 


0   i=l 


£  0,    6x,(T)  =  6S, 


i    1 


(C-16) 
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where  6S  represents  the  variation  In  the  objective  function. 

From  the  Taylor  Series  expansion  of  a  function  about  a  point 
we  have, 

f(xrV  ••"  xn>  =  f(al'R2 V 

n  i  . 

+  higher  order  terms 

where  the  second  order  partial  derivatives  are  assumed  to  exist. 
Using  this  result,  the  factor  within  brackets  in  the  second 
term  on  the  right  hand  side  of  equation  (C-l?)  becomes 

f  (x  +  6x;  9  +  68)  -  f  (xj  6) 

n  xf   (x,  G  +66) 

=  f,(xj  e  +40)  -  f.(xj  e)  +  s  2_i         6x 

+  higher  order  terms      (C-19) 

Recalling  equation  (C-l?)  we  have 

In  fe.  T  n 

as  =    /    z    -~Hx  at  +    /    j;  z  [f  (x+6xj  ?+6e)  -  f,(xje)]dt 

0      1=1     dL       1  0      1=1   *      1  1 


and  finally  utilizing   equations    (C-9)    and   (C-19)   gives 

T  fi     n         Jf  T  n 

6S  ,-    -   /      5.      i;   z     r^6xdt  +      /      £   z,[f,  (x,   +69 )    -   f,(x;"e) 
0     1=1  J=l  J  ^    i     1  0     1=]    *     i  1 


n    }f 


+     £     \-~(x!    8   +   68)6x  Idt 
J=l  «xj  J" 

T     n 
/        I.      2  Tf    (X!   ?  +   68)    -   f,  (xj    e")ldt 
0        1=1      1      *•  l 
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T  n  n    if, if  _  _ 

.;-   /  £  £  Z,  £-i(xj8+6e)  -  2_i.(xj6)  6x  dt  (C-20) 

0   1=13=1  *•  dxj  d'j        J 

At  this  juncture  we  shall  impose  a  restriction  on  the  functions 

f  whereby  we  shall  only  consider  functions  that  are  linear  in  x 
1 

and  that  contain  6  in  a  separable  fashion;  i.e.  the  performance 

equations  shall  be  of  the  form  shown  in  equation  (C-3). 

n 
£ 
1=1   1J  1 


f,  =  £  a   x  +  Tfte)  (c-21) 

_n   iJ   J   ' 


where  a   is  not  necessarily  constant. 

For  this  class  of  functions  we  obviously  have 

^-  fjTx;  tf  +  66)  =  JL  f^xj  V)  (C-22) 

J  J 


and  this  fact  causes  the  second  integral  in  equation  (C-20)  to 

vanish j  equation  (C-20)  then  reduces  to 

T  n      _  _ 
/   £  B.CM*i  6  +  66)  -  f  (xi  6)]dt 
3    U1   1 


63  = 

0 


T 
/  [H(x;  G  •(•  66)  -  H(x;  6)]dt  (C-23) 

0 


Since  S  is  to  be  maximized,  it  is  necessary  that  6S  be  zero  along 
the  optimal  trajectory  for  all  free  variations  6  6  and  that  6S 
be  negative  for  variations  66  at  the  boundary  of  the  constraints 
as  expressed  by  equation  (C-5).   Therefore 

63  <  0  (C-2^) 
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and  hence  by  virtue  of  equation  (0-23)  it  follows  that 

6H  =  H(xj  "5  +  6  6)  .  H(X|  I*)  <  0  (C-25) 

in  the  total  interval 

0  <  t  <  T 

Also,  since  equation  (C-25)  holds  for  any  small  perturbation  6  8, 
it  is  necessary  for  the  Kamiltonian  to  take  on  a  maximum  value 
when  S  is  to  be  maximized  and  on  the  other  hand  it  is  necessary 
for  the  Hamiltonian  to  be  a  minimum  when  S  is  to  be  minimized. 
An  interesting  geometrical  interpretation  of  the  Maximum 
Principle  may  be  obtained  if  one  considers  the  fact  that  H  is 
really  the  scalar  product  of  the  z(t)  and  f(x,  6)  vector  functions 
and  that  the  6(t)  vector  should  be  so  mainpulated  so  as  to 
maintain  this  scalar  product  at  a  maximum  or  minimum  according 
as  we  wish  to  maximize  or  minimize  the  objective  function. 

C-2.2.   FINAL  TIME  FIXED  AND  STATE  VARIABLES  WITH  EQUALITY 
CONSTRAINTS 

We  shall  now  consider  the  case  where  the  final  state 

variables  have  equality  constraints  as  follows, 

F  [x(T)]  =  0        a  =  1,  2 p  (C-26) 

a  P  <   n 

Here  we  shall  employ  the  variational  technique  using  Lagrange 

Multipliers.  Consider  the  objective  function  of  the  form 

n  p 

S=  T.     C   x(T)  .;•   I:  \_   F  O(T)]  (C-2?) 

1=1  *  *      a=l     a 


1?0 


instead  of  the  original  form 

n 
S  =   X  C   x.  (T) 
1=1   1   1 


and  v:e  desire  to  maximize  S. 

We  see  from  equations  (C-l?)  through  (C-20)  that 

n       ,T     T  n 

£  Ux.    =   /  T.   z,[f.(x;  "0  +  66)  -  f  (xi  6)jdt 
1=1  i   l  lo   0  1=1  x      x  1 

Tn  n   Nf  (xi'e+ee)  Kitx'si") 
0  1=13=1  *  oxj        6  j 


dt 


(C-28) 

In  order  that  the  left  hand  side  of  equation  (C-28)  represent 
the  variation  in  the  objective  function,  the  boundary  conditions 
of  the  z(t)  vector  have  to  be  determined  so  that  equation 
(C-25)  holds. 

From  equation  (C-27)  we  have 


n  P       n  £  P0. 

6S  =      L  C    6x   (T)   +     £        EX 
1=1  *      x  o= 


2  X     - — ~6x,(T)   +   higher        (C-29) 

1U°«V    '  order  terms 


If  the  constraints  F  have  finite  values  for  second  order 
derivatives,  then  in  the  neighborhood  of  optimal  values  we  may 
choose  X  _  such  that 

B        p     ^F       I 

6S  =  t   [C  -i-  5J   X  ^-aj  6X  <C-30) 

1=1   i   o-l   °'  ^    1  I  t  =  T 

Now  comparing  the  left  hand  side  of  equation  (C-28)  to  equation 
(C-30).  it  Is  easy  to  see  that 


P     >F 

x  SLai 


z  (t)  =  [c,  +    x    x   :L-2n 

l        i   0=1  *  &Vt,.T 


(0-31) 
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Therefore,  we  now  have  n  equations  from  equation  (C-l),  n 
equations  from  equation  (C~31)  together  with  p  equations  from 
equation  (C-26)  with  the  help  of  which  we  can  solve  for  the  2n 
unknown  components  of  x(t)  and  z(t)  and  the  p  unknown  components 
of  Xa,   Equation  (C-31)  Is  equivalent  to  the  transversality  con- 
dition of  the  final  values  in  the  calculus  of  variations.   A 
special  case  of  the  above  problem  is  when  only  some  of  the  final 
values  of  the  state  variable  are  fixed,  say 

x  (T)  =  x^     i  =  1,  2,  ....  m 

and  the  rest  of  the  final  values  are  still  free.   In  this  case, 
equation  (C-20)  becomes, 

F  [x(T)]  =  x ,  (T)  -  X?  =»  0  (C-32) 

o.        1       i 

Substituting  equation  (C-32)  into  equation  (C-31)  yields 

2  (T)  m  lCx   +  Xt]        1  «  1,  2,  ....  B  (C-33) 

z AT)   =  Ct  I  =  m+1,  m+2,  ,..,  n 

Next  let  us  consider  the  case  where  the  constraints  are  imposed 
on  the  initial  values  of  the  state  variables  In  the  form 

F  [x(0)]  *  6  p  «  1,  2 p  (C-3'l) 

P  <;  n 

From  equation  (C-28)  it  follows  that  at  the  optimal  condition 


1=1  *   1  lt-0 
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If  the  initial  values  arc  preasslgned,  then  fix  (0)  =  0  holds, 
and  equation  (C-18)  is  thereby  satisfied!  But  if  the  initial 
values  have  a  restriction  as  in  the  form  of  equation  (C-3'i-). 
then  in  general 


z,    fix 


1=1 


n  n 

£  C      fix  ,  (T)    -     £   z,    fix. 

1=1    1  X  lal    ■ 


i|t=0 


63 "  Ji  zi £x^ 


Since  at  the  optimal  condition  we  have  6S  =  0  it  follows  that 


5Jz  fix,     =0  (C-36) 

1    A  |t=0 

It  should  be  noticed  that  the  variation  fix  (0)  in  equation 

(C-36)  is  under  the  restrictions  of  equation  (C-3'O.   Since  small 

perturbations  around  optimal  conditions  have  been  considered,  the 

variations  fix  (0)  lie  on  the  pianos  tangent  to  the  surfaces  of 

equation  (C-34).  Equation  (C-36)  is  the  transversal ity  condition 

of  the  Initial  values. 

C-2.3.   FINAL  TIME  UNSPECIFIED 

The  previous  cases  considered  were  characterized  by  a 
fixed  final  time  T,  we  shall  now  consider  the  case  where  the 
final  time  T  is  free  and  is  to  be  determined  along  with  the 
optimal  state  vector  x(t)  and  the  optimal  decision  vector  ¥(t). 
This  is  what  is  known  as  the  time  optimal  problem. 

If  the  final  time  T  appears  In  the  objective  function,  we 


introduce  an  n+1  state  variable  x    such  that 

n+l 
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n±i=  i,       x   (o)  =  0  (C-37) 

at  n*l 


and  x   (T)  is  identical  to  T  Itself, 
n+l 


n+l  >, 

6S  = 


1=1  &V  '    i 


n+l       P      X  F„   dx, 
i=l  i   a-l  a  ixi(T)  ^ 


t=T 


Then  it  follows  from  the  above  equation  and  equation  (C-30)  that 

n+l       p 
6S  =  £  [C  +  T.     \ 
1=1   3   a  " 

n+l    dx, 

i=l  i  dt     U+T 

Since  we  want  to  maximize  S,  6S  <  0.   But  on  the  other  hand  6t 

may  take  positive  or  negative  values  and  therefore  as  S  attains 

its  maximum  value  63  must  equal  zero.   It  therefore,  follows 

that 

n+l    dx.j 
H(T)  «=  L  z,  -rp     -_-  0  (C-39) 

Also  since  H  has  a  constant  value  throughout  the  transient  state, 
it  follows  that 

H  =.-.  0,       0  <  t  <  T 

is  a  necessary  condition  for  optimality  when  the  final  time 
T  is  unspecified, 


l?k 


ZA     EXISTENCE  THEOREMS   C?  MAEKUS   AND  LEE  [13] 
Theorem  1.     Given  the  control  problem 

(a)  x    =  fjL(t,   Xj,   x?,    ,,.,    xn,   Uj ura) 

b  8j(ti   1)  +  hj(ti   x)  U,  1  *  1,   2,    ...,   n 

j   =   1 ,    2 ,    .  .  .  ,    in 

Kith  « 

1  dM*-  >.)    Xh,(tr  x) 

g   (t,    x),    hAJ(t,  x), i -,r-i ,      k  =  1,    2 n 

i  1  hk  }xk 

continuous  in  all  their  variables, 

(b)  a  nonempty,  convex,  compact  restraint  set  X),  in  R  , 

(c)  the  initial  point  x  in  R  , 

(d)  the  continuously  moving  nonempty  compact  target  set  G(t) 
on  a  finite  interval  yo  <  t  <  v  , 

(e)  the  cost  functional 

h 

C(u)  b   /   f  (t,  x(t),  u(t))dt 

t      ° 
o 

where  f-(t,  x,  u)  ~   g  (t,  x)  +  h  (t,  x)u  and  the  functions 

g  (t,  x)  and  h  (t,  x),  J  =  1,  2,  ..,,  m  are  continuous  In  all 
o  o 

their  variables. 

Assume  the  set  A  of  controls  with  responses  traveling  from 

x  to  G  is  such  that 
o 

(A)  A  Is  non  empty 

(B)  there  exists  a  real  bound  B  <  »  for  all  responses  x(t) 

.  n 

corresponding  to  u(t)  in  A ,  that  is,  |x(t) I-  E  (1  (t)j 

1  i-l'  i   ' 
<  B  uniformly  for  all  responses. 


i?; 


Then  there  exists  an  optimal  control  u  (t)  in  A. 

Remark i   Hypothesis  B  of  theorem  1  la  satisfied  for  x  In  R  and 

u  in  IT. If 


or  If 


f  ( t ,  x ,  u )  I  <  a ,      1  =  1 ,  2 ,  ....  n 
^(t,  x,   u)j 


i*J 


<  a, 


1,  2, 


for  some  real  a.  Hypothesis  A  of  theorem  1  deals  with  the 

question  of  whether  or  not  the  system  can  be  moved  from  x     to 

G(t)  in  a  finite  time.   The  set  of  initial  points  C   from  which 

o 

the  system  can  be  moved  to  G(t)  in  a  finite  time  Is  called  the 
domain  of  controllability  and  is  the  subject  of  the  next  theorem. 
Theorem  2.   Consider 


dx. 

dt  -   Vxl' 


K  ,  ix,  ,  •  •  •  i  u  ) , 

n'  T.'       m" 


1  =  1 n 


■where  f(x,  u)  and  A£(x,  U),  and  -«J-(x,  u)  are  continuous  in  R  x 
SI.      The  control  restraint 51 Cft  contains  the  origin  in  its 
interior. 
Assume I 

(1)   f(0,  0)  =  0 

{?.)      There  exists  a  vector  v  £  Rm  such  that  Bv,  ABv,  ...,  A   Bv, 
are  linearly  independent,  where  A  ^  f— i.(0,  0}  and  E  =  r~-=-(0,  0), 


*u 


1.  2, 


, ,  n  j  k  a  1 ,  2 , 


k 


Then  there  exists  a  neighborhood  U  c  R  of  the  origin  such  that 


raoh  point  x  c  U  can  be  steered  to  the  origin  in  R  in  a 
0 


1?6 


finite  time  interval,  using  a  measurable  control  function  u(t) 
viith  graph  inJ^.  Remarki   Hypothesis  (?)  is  satisfied  if 

det  [Br,  ABv,  ...,  fn~  Bv]  J   0« 
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ABSTRACT 

Host  of  the  literature  dealing  with  control  of  distillation 
columns  considers  simplified  models  of  complete  mixing  to 
represent  the  liquid  phase  mixing  on  the  distillation  trays. 
In  this  work  a  two  tray  distillation  column  is  considered  wherein 
the  top  tray  is  described  by  the  so-called  mixing  pool  model 
consisting  of  two  tanks  In  series.  Also  for  the  sake  of 
comparison  a  conventional  distillation  column  with  complete 
mixing  on  the  top  tray  Is  also  considered.  The  control  problem 
is  defined  thusi  the  system  suffers  a  disturbance  through  the 
feed  composition  which  in  turn  displaces  the  overhead  distillate 
composition  from  its  steady  value.   Determine  a  control  policy 
for  the  overhead  reflux  rate  such  that  a)  the  overhead  distillate 
composition  will  be  returned  to  its  steady  state  value  in  the 
shortest  possible  time  b)  the  overhead  distillate  composition 
will  be  returned  to  its  steady  state  value  and  in  doing  so  Its 
deviation  from  the  steady  state  value  will  be  minimum  in  a 
least  squares  sense.   Both  problems  a)  and  b)  are  applied  to 
both  systems  described  above  and  their  corresponding  control 
policies  are  obtained.   It  is  found  that  the  change  in  the  model 
has  a  significant  effect  on  the  control  policy  and  on  the  methods 
of  obtaining  the  control  policies.   Various  extensions  to  the 
present  problem  are  also  suggested. 


